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ABSTRACT 
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The  present  work  is  a  Kalman  filter  study,  in  indirect 
feedback  configuration,  for  a  proposed  integrated  inexpen¬ 
sive  Inertial  Navigation  System/Global  Positioning  System 
(I.N.S./G.P.S. ) . 

A  one  nautical  mile  per  hour  ,  local-level  ,  two- 
accelerometer  I.N.S.  is  used  where  the  errors  are 
represented  by  a  7  state  linear  model . 

G.P.S.  is  assumed  to  provide  four  range  measurements 
from  an  equal  number  of  satellites  with  the  best  relative 
position  among  those  in  view. 

I.N.S.  error  analysis  showed  error  dependence  on 
Schuler  frequency  and  that  it  was  possible  to  neglect 
Foucault  modulation  for  navigation  purposes. 

The  present  I.N.S./G.P.S.  system  has  been  shown  to  be 
quite  effective  since  the  navigation  errors  are  reduced 
quickly  for  both  short  and  long  term  periods  without  any 
divergence. 
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I.  INTRODUCTION 


A.  OVERVIEW  OF  AN  INERTIAL  NAVIGATION  SYSTEM  (I.N.S.) 

A  conventional  gimballed  inertial  measurement  unit 
consists  of  a  platform  suspended  by  a  gimbal  structure  that 
allows  three  degrees  of  rotational  freedom  [Ref.  1,2,7]. 

The  outermost  gimbal  can  be  attached  to  the  body  of  some 
vehicle  and  allow  that  vehicle  to  undergo  any  change  in 
angular  orientation  while  maintaining  the  platform  fixed 
with  respect  to  some  desired  coordinate  frame. 

Gyros  mounted  on  the  platform  sense  the  angular  rate  of 
the  platform  with  respect  to  inertial  space  and  their 
outputs  are  sent  through  electronics  to  tne  torque  motors  on 
the  gimballed  structure,  commanding  them  to  maintain  a 
desired  platform  orientation  regardless  of  the  orientation 
of  the  outermost  gimbal  which  remains  fixed  to  the  body. 

Feedback  control  loops  that  keep  the  gyro  outputs 
nulled,  will  maintain  at  the  same  time  the  platform  fixed 
with  respect  to  the  inertial  space.  These  feedback  loops 
are  such  that,  in  practice,  the  platform  orientation  is  kept 
essentially  stable  regardless  of  the  most  violent  vehicle 
maneuvering.  Additional  (computed)  inputs  can  be  added  to 
the  above  feedback  loops  to  maintain  some  other  orientation, 
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such  as  North-East-Down,  corresponding  to  the  current 
location  of  the  vehicle. 

Accelerometers  mounted  on  the  platform  can  provide  the 
vehicle's  acceleration  with  respect  to  the  known  set  of 
reference  coordinates.  in  fact,  specific  force  is  measured 
by  the  accelerometers  so  that  local  gravity  must  be  computed 
and  appropriately  subtracted  from  these  sensor  outputs 
in  order  to  obtain  a  measurement  of  actual  vehicle 
acceleration  . 

The  vehicle's  velocity  and  position  are  obtained  by 
integration  of  the  above  acceleration  measurement  signals. 
Attitude  information  as  well  as  translational  information  is 
provided  by  the  I.N.S..  A  typical  gimballed  inertial 
measurement  unit  [Ref.  2]  is  shown  in  Fig.  1. 

B.  OVERVIEW  OF  THE  G.P.S. 

The  Global  Positioning  System  (GPS)  is  a  satellite 
navigation  system  currently  under  development.  It  will 
consist,  according  to  today's  available  information,  of  18 
satellites  placed  in  groups  of  six  in  each  of  three 
different  circular,  12  hour  orbits  at  an  altitude  of  10,900 
N.M.  inclined  63°  to  the  equator  and  spaced  120°  apart. 

The  satellites  will  broadcast  pseudo-random  noise  codes 
(codes  P  and  C/A)  and  ephemerides  on  two  L-band  signals  to 
users  worldwide  in  such  a  way  that  each  satellite  signal  can 
be  distinguished  from  the  others  by  the  user.  A  user  will 
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Figure  1.  Typical  Gimbaled  Inertial  Measurement  Unit  [Ref.  2] 


be  equipped  with  a  small  receiver  (G.P.S.  user  equipment) 
which  measures  the  pseudo-range  and  pseudo-range  rate  from 
the  user  to  the  satellite. 

By  means  of  a  correlator-detector  the  time  (phase)  shift 
between  each  satellite  signal  and  the  user's  unsynchronized 
clock  will  be  measured  in  his  receiver  to  provide  an 
indication  of  the  range  from  the  satellite  to  the  user. 
Typically,  four  satellite  signals  may  be  received  simul¬ 
taneously  by  the  user  equipment. 

The  phases  of  the  NA  VSTAR/G .  P.  S.  system  are  snown  in 
Fig.  2  [Ref.  8]. 

C.  I.N.S.  OPTIMAL  AIDING 

Once  we  have  available  a  typical  inertial  measurement 
unit,  or  the  inerti3i  navigation  system  as  a  whole,  the 
question  naturally  arises:  why  does  this  system  requires 
optimal  aiding  by  other  navigational  sensors?  The  answer  to 
thi3  question  is  given  in  Ref.  2  and  here  we  present  the 
concepts  only. 

Due  to  the  tight  control  loops  supporting  the  I.N.S. 
very  good  high  frequency  information  is  provided.  However, 
because  of  gyro  characteristics,  the  system  drifts  at  a  slow 
rate  so  that  the  long  term  (low  frequency  content)  of  the 
data  is  poor.  It  is  well  known  that  all  inertial  systems 
have  position  errors  that  grow  slowly  with  time  and  these 
errors  are  unbounded. 


POSITIONING  SYSTEM  SCHEDULE 


NAVSTAR/G.P.S.  Phases 


As  opposed  to  an  I.N.S.  which  can  be  classed  as  a  "one 
nautical  mile  per  hour  system"  due  to  the  associated 
position  error,  most  other  navigation  aids  provide  very  good 
low  frequency  information  but  subject  to  considerable  high 
frequency  noise,  due  to  instrument  noise,  atmospheric 
effects,  antenna  oscillation,  unlevel  ground  effects  and  so 
for  th . 

One  would  want  to  combine  the  available  information  from 
an  I.N.S.  and  other  external  sources  in  an  optimal  manner  if 
possible  so  that  one  can  obtain  efficient  estimates  of 
navigation  parameters  that  are  best  with  respect  to  some 
well  defined  criterion.  Such  an  optimal  approach  is 
provided  by  the  Kalman  filter  approach  which  is  briefly 
discussed  next. 

D.  KALMAN  FILTER 

The  Kalman  filter  is  an  optimal  recursive  data  process¬ 
ing  algorithm  located  in  the  on-board  computer  or  central 
processor  that  uses  sampled  data  with  sample  period  on  the 
order  of  5-60  seconds,  to  maintain  estimates  of  approxi¬ 
mately  60-70  state  variables.  The  filter  combines  all 
available  measurement  data  with  prior  knowledge  of  the 
system  and  measuring  devices  to  produce  an  estimate  of  the 
system  states  in  such  a  manner  as  to  statistically  minimize 
the  resulting  errors.  In  more  easily  understood  terms  the 
filter,  or  computer  program,  uses  the  statistical 
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characteristics  of  the  errors  in  both  the  inertial  naviga¬ 
tion  components  and  the  external  information  providing  the 
best  estimate  possible,  subject  to  certain  modeling 
assumptions . 

The  filter  will  act  to  optimize  the  attitude,  position, 
and  velocity  information  accuracy  by  weighting  each  data 
source  heavily  in  the  frequency  ranges  where  it  provides 
more  accurate  information,  and  suppressing  it  in  the  region 
wnere  it  ts  less  accurate.  The  inertial  system  provides 
good  high  frequency  information  but  it  drifts  slowly  and 
therefore  exhibits  poor  low  frequency  performance.  On  the 
other  hand,  the  external  aids  (such  as  G.P.S.)  generally 
exhibit  good  low  frequency  information  but  are  subject  to 
high  frequency  noise.  Therefore,  the  filter  will  use  the 
good  low  frequency  external  (G.P.S.)  information  to  damp  out 
the  slowly  growing  errors  in  the  inertial  system. 

1 .  Type  of  Filter  Implementation 

There  are  two  very  important  aspects  of  implemen¬ 
tation  of  a  Kalman  filter  in  conjunction  with  an  inertial 
system  [Ref.  2]. 

a'  Total  state  space  (direct)  versus  error  state 
space  (indirect)  formulation,  and 

b)  Feedforward  versus  feedback  mechanizations. 

In  the  indirect  formulation  the  errors  in  the  I.N.S. 
indicated  position  and  velocity  are  among  the  estimated 
variables  and  each  measurement  presented  to  the  filter  is 
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the  difference  between  the  I.N.S.  and  the  external  source 
(G.P.S.)  data.  The  I.N.S.  itself  follows  the  high  frequency 
motions  of  the  vehicle  very  accurately,  and  there  is  no  need 
to  model  these  dynamics  explicitly  in  the  filter  but  the 
dynamics  upon  which  the  filter  is  based  is  a  set  of  inertial 
system  error  propagation  equations,  which  are  relatively 
well  developed,  well  behaved,  low  frequency,  and  very 
adequately  represented  as  linear  [Ref.  2,  pp .  296]. 

The  indirect  feedback  conf iguration  is  considered 
where  the  Kalman  filter  generates  the  estimates  of  the 
errors  of  the  I.N.S,.  and  feeds  bacx  these  errors  to  the 
I.N.S.  to  correct  it.  By  this  configuration  we  use  the  two 
major  advantages.  First  that  the  I.N.S.  errors  are  not 
allowed  to  grow  unchecked  and  the  adequacy  of  a  linear  model 
is  enhanced.  Second  is  the  fact  that  many  of  tne  predicted 
error  states  which  at  next  time  sample  time  are  zero  ,  need 
not  be  computed  explicitly. 

The  indirect  feedback  configuration  of  the  Kalman 
filter  is  shown  in  Figure  3.  In  Ref.  2  there  is  explicitly 
documented  the  discussion  of  the  Kalman  filter  configura¬ 
tions  and  mechanizations  with  their  advantages  and  more 
comments . 

2.  Assumptions 

The  Kalman  filter  can  be  shown  to  be  the  best  filter 
of  any  possible  form  based  on  the  following  three 
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assumptions:  a  linear  system,  white  noise  drivers,  and  the 
Gaussian  distribution  of  noise. 

Although  the  system  itself  may  be  nonlinear, 
formulation  of  an  approximate  linear  error  state  space  model 
makes  linear  analysis  possible.  The  justification  for  the 
linear  model  is  based  on  two  points.  For  the  aided  I.N.S. 
case  the  use  of  linear  error  state  space  models  yields  a 
very  adequate  representation .  The  techniques  of  linear 
system  analysis  are  also  well  developed  and  better 
understood  than  those  of  nonlinear  analysis. 

The  white  noise  assumption  implies  that  the  noise 
is  not  correlated  in  time  and  also  has  equal  power  at  all 
frequencies.  If,  in  fact,  a  time  correlated  noise  is 
required  to  adequately  model  the  system,  it  can  be  produced 
by  passing  white  noise  through  a  linear  shaping  filter.  The 
system  can  then  be  modeled  with  an  augmented  state  variable 
as  a  linear  system  driven  by  white  noise. 

Gaussianess  pertains  to  the  distribution  of 
amplitudes  of  the  noise  and  implies  that  at  any  single  point 
in  time  the  probability  density  function  of  the  amplitude 
takes  on  the  shape  of  a  normal  bell-shaped  curve.  The 
assumption  of  Gaussian  noise  amplitude  is  justified  by  the 
fact  that  the  system  or  measurement  noise  is  typically 
caused  by  a  number  of  sources.  It  can  be  shown 
mathematically  that  when  a  number  of  independent  random 
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variables  are  added  the  resultant  effect  is  very  nearly  a 
Gaussian  probability  density  even  though  the  individual 
densities  are  not  Gaussian  [Ref.  93* 

Under  the  above  mentioned  assumptions  of  a  purely 
white  Gaussian  noise,  the  first  two  moments  specify  the 
entire  shape  of  the  density  describing  the  noise,  and  the 
mathematics  of  the  problem  are  greatly  simplified. 

In  Appendix  A,  a  simple  example  of  a  Kalman  filter 
application  to  a  radar  position-aided  I.N.S.  is  given  in 
order  to  make  easier  the  understanding  of  Kalman  filter 
operation . 

Finally,  information  about  the  G.P.S.  satellites' 
geometry  and  their  observability  is  given  in  Appendix  B. 
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II.  KALMAN  FILTER  EQUATIONS 


A.  GENERAL 

The  design  of  a  Kalman  filter  and  especially  the 
integrated  I.N.S.  Kalman  filter  design  requires  extensive 
computer  simulation.  This  chapter  of  the  work  is  a 
presentation  of  the  equations  which  are  required  not  only 
for  the  mechan i za t ion  of  the  filter  but  also  tnose  which  are 
necessary  to  simulate  the  dynamics  of  any  user  (aircraft, 
missile)  . 

The  principal  tool  used  for  the  solution  of  this 
specific  and  other  similar  problems  is  the  very  common 
method  of  aov  ar iance  analysis.  It  is  known  that  the 
covariance  is  a  measure  of  the  uncertainty  in  v.ne  knowledge 
of  the  true  values  of  the  state  vector  components.  In  this 
work,  as  the  covariance  matrix  is  concerned,  the  off 
diagonal  terms  are  assumed  to  be  zero  initially  and  initial 
conditions  on  the  diagonal  elements  are  arbitrarily  taken. 
The  covariance  matrix  of  both  the  system  and  the  filter  are 
propagated  forward  in  time  by  numerical  integration 
techniques . 

The  adjustment  of  the  values  of  the  state  variables,  to 
those  of  the  best  estimate  obtained  with  the  Kalman  filter, 
is  achieved  when  a  control  is  applied  to  the  system  after 
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the  specified  update  time  is  reached  and  the  best  estimates 
of  the  states  have  been  determined.  The  square  root  of  the 
individual  diagonal  elements  of  the  system  covariance  matrix 
(RMS  values)  are  plotted  as  a  function  of  time  to  provide 
the  performance  of  the  filter.  For  this  study  the  plots  are 
also  utilized  to  determine  the  error  contribution  associated 
with  each  modeled  error  source.  Furthermore,  the  error 
statistics  are  propagated  which  means  that  the  standard 
deviation  of  the  noise  value  is  supplied  whenever  a  noise  is 
required . 

One  attribute  of  covariance  analysis  is,  that  under  the 
assumptions  stated  in  Chapter  One,  i.e.,  the  linearity  and 
white  Gaussian  noise,  the  covariance  is  independent  of  the 
actual  measurement  values  and  can  be  computed  through 
generating  a  sample  sequence  of  measurements.  And  as  a 
matter  of  fact  this  method  is  easier  to  handle  and  work  with 
than  the  correspond ing  Monte  Carlo  type  simulation. 

B.  SYSTEM  MODEL  EQUATIONS 

The  differential  equations  that  describe  how  the 
inertial  navigator  errors  propagate  with  time  are  the  basic 
equations  used  in  this  process .  These  equations  are 
formulated  in  a  set  of  first  order,  linear  differential 
equations,  driven  by  white  Gaussian  noise  for  the  reasons 
described  previously  in  this  work. 
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Linear  measurements  corrupted  also  by  white  Gaussian 
noise  are  made  upon  the  actual  system  variables.  It  is 
furthermore  assumed  that  the  equations  which  represent  a 
detailed  model  of  the  system  are  of  the  form: 
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characteristics  of  zero  mean  and  variance: 


T  Q_(i>  for  i  s  J 

EC w( i) w( j) 1 1  =  (2  ) 

0  for  i  i  j 


where  the  indices  i  and  j  are  instants  in  time. 


The  observations  which  are  obtained  from  external 
references  and  in  our  case  of  study  from  the  G.P.S.  can  be 
described  by  the  following  linear  measurement  vector 
equation : 
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where : 

z_  is  a  q  vector  of  measurements 

““  3 

Hs  is  a  q  x  n1  matrix  of  measurements 

v3  is  a  q  vector  of  white  noise  inputs  with  the 
characteristics  of  zero  mean  and  variance: 

P  H  ( i)  for  i  =  j 

E[vU)v(j)  j  =  3  (4) 

0  for  i  i  j 

A  further  assumption  for  the  study  is  that  the  system  noise 
w  and  the  measurement  noise  v  are  uncorrelated  for  all  time, 
i  .e  .  , 

E[w(i)v(j)T]  =  0  for  all  i,j  (5) 

C.  FILTER  EQUATIONS 

The  equations  discussed  above  are  assumed  to  be  a 
complete  and  accurate  mathematical  description  of  the  G.P.S. 
aided  inertial  navigation  system  dynamics  and  measurement 
equations  for  the  purpose  of  simulation.  They  also 
constitute  a  set  of  equations  which  would  be  utilized  in  the 
design  of  a  fully  optimal  Kalman  filter. 

In  our  case  of  study  as  also  in  general  a  suboptimal  or 
reduced  order  filter  design  is  obtained  by  reducing  the 
dimension  of  the  state  vector  due  to  the  computational 
burden  of  the  fully  optimal  filter.  The  states  that  are 
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eliminated  are  those  that  affect  the  accuracy  the  least  of 
the  mathematical  description  of  the  aided-I.N.S.  The 
designed  suboptimal  filter  can  be  implemented  with  the  on¬ 
board  computer  (aircraft  or  missile). 

The  equations  which  describe  the  suboptimal  filter  are 
of  the  form: 

if  =  Ff*f  +  Gf«f 

where 
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The  equation  for  the  filter  measurement  is : 

If  s  HjOCf  v^.  (7) 

where : 

2^  is  a  q  vector 
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p 


is  a  q  x  q  matrix  of  measurements 

v^  is  a  q  vector  of  white  noise  inputs  with  the 
characteristics  of  zero  mean  and  variance: 

T  Rf(i)  for  i  =  j 

E[v-(i)v-(  j)  3  =  (3) 

0  for  i  i  j 

The  filter  propagation  and  update  equations  based  on  tne 
above  models  are  then  given  below. 

At  measurement  times  (update): 

Kf  -  ff'HfhVfV  *  afr' 

P  p  —  P  f  —  K  f  d  P 
_xf+  =  +  Xf[zs  " 

and  between  measurements  (extrapolate): 

If  =  Fflf 

Pf  =  FfPf  +  PfFfT  ♦  GfQfGfT 

where : 

is  an  ng  vector  denoting  the  best  estimate 
Pj*  is  the  covariance  matrix  of  the  filter 
Kf  is  the  matrix  of  Kalman  gains 


(9) 

(10) 

(11) 

(12) 

(13) 
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:  is  a  q  vector  of  the  actual  values  of  the 

— s 

measurements  taken 

+  superscript  indicates  the  time  instant  just  after 
update 

-  superscript  indicates  the  time  instant  just  prior 
to  update 

T  superscript  denotes  the  transpose  matrix  or  vector 
superscripted  . 

The  filter  subtracts  from  the  actual  taken  measurement 
_zs  the  best  prediction  of  its  value  before  the  actual 
measurement  is  taken,  i.e.,  the  value  of  This 

difference  is  then  passed  through  an  optimal  weighting 
matrix  and  used  to  correct  x^”,  the  best  prediction  of 
the  state  at  the  time  instant  before  the  measurement  is 
taken.  This  process  gives  the  best  estimate  after  update. 
This  estimate  is  propagated  to  the  time  of  the  next 
measurement  sample  according  to  equations  (12)  and  (13). 

The  above  recursive  relationships  are  solved  based  on 
initial  conditions  of  an  assumed  Gaussian  density  which 
describe  the  a-priori  knowledge  of  the  state  as: 

2(0)  =  xQ  (14) 

and  P(0)  =  PQ 
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The  Kalman  filter  conditioned  on  the  actual  measurements 
taken,  propagates  the  conditioned  probability  density  of  the 
desired  states.  The  probability  density  function  of  a 
Gaussian  noise  amplitude  takes  on  the  shape  of  a  normal 
bell-shaped  curve.  The  assumption  of  Gaussian  noise 
amplitude  is  well  justified  by  the  fact  that  a  system  or 
measurement  noise  is  typically  caused  by  a  number  of  small 
sources  and  according  to  the  Central  Limit  Theorem  it  can  be 
shown  mathematically  that  when  a  number  of  random  variables 
are  added  together,  the  summation  is  a  random  variable  whose 
density  is  nearly  a  Gaussian  probability  density,  regardless 
of  the  shape  of  the  densities  of  the  individual  random 
variables.  Furthermore,  the  use  of  Gaussian  densities  makes 
the  mathematics  easier  to  handle  and  tractable.  It  is  known 
that  a  Gaussian  density  is  completely  determined  by  its 
first  and  second  order  statistics,  i  .e .  ,  the  mean  and  the 
variance.  Thus,  the  Kalman  filter,  which  propagates  the 
first  and  second  order  statistics,  includes  all  information 
contained  in  the  conditional  probability  density  mentioned 
above  [Ref.  2:  pp.  3-9]* 

The  mean  of  a  density  function  or  its  expectation  u,  is 


defined  as 


and  it  is  interpreted  as  the  weighted  average  of  the  values 
of  x,  using  the  probability  density  function  f(x)  as  the 
weighted  function .  All  the  Gaussian  white  noise  inputs  in 
this  study  are  assumed  to  have  zero  mean. 

The  variance  of  a  density  function tor  the  square  of  the 
standard  deviation  c,  is  defined  as: 

V 3 r [ x ]  =  a2  =  f  ( x  -  y)2f(x)dx  (17) 

•  CO 

and  it  is  interpreted  as  the  weighted  average  of  the  values 
of  ( x  -  u)~;  thus,  a  is  a  measure  of  the  density  spread  and 
a  direct,  measure  of  the  uncertainty  since  the  larger  j  is, 
the  broader  the  probability  pea'x  is,  spreading  the  proba¬ 
bility  weight  over  a  larger  range  of  x  values.  For  the 
example  of  Gaussian  density,  65.3%  of  the  probability  weight 
is  contained  within  the  band  of  a  units  to  each  side  of  the 
mean  u,  which  represents  the  area  under  the  normal  bell¬ 
shaped  curve  between  tht  values  of  -a  and  +  a  and  95.4%  of 
the  probability  weight  is  contained  between  the  values  of 
-2  a  and  +2a.  Since  in  this  study  vectors  are  used  instead 
of  3calars,  the  above  equations  (16)  and  (17)  which  give  the 
first  and  second  order  statistics  respectively  for  the 
scalar  case,  must  be  extended  for  the  vector  case  as 
follows  : 
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III.  I.N.S.  ERROR  ANALYSIS 

A.  GENERAL 

One  general  approach  to  determine  the  navigation  error 
caused  by  individual  sources  of  error  is  to  simulate  the 
inertial- navigation -system  nonlinear  equations  an j  sources 
of  error  and  compare  the  navigation  outputs  with  the 
simulated  true  position--the  differ. ..ce  being  the  navigation 
error.  However,  this  is  not  the  approau-n  used  here. 

Assuming  that  tne  position  errors  are  small  compared  with 
earth  radius,  that  the  velocity  errors  are  small  compared 
with  orbital  velocity,  and  that  the  alignment  errors  are 
small  compared  with  1  radian  (or  3437  arc-min)  it  can  be 
demonstrated  tnat  the  propagation  of  errors  in  an  inertial 
navigation  system  is  very  accurately  governed  by  a  set  of 
corresponding  linear  differential  equations.  Therefore, 
most  inertial-navigation-system  error  analyses  are  conducted 
working  directly  with  a  set  of  linear  error  differential 
equations . 

Sets  of  error  equations  have  been  developed  for  various 
I.N.S.  conf igurations  in  the  context  of  a  particular 
application.  As  a  consequence,  many  sets  of  I.N.S.  error 
equations  have  been  developed  for  the  various  broad  classes 
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of  mechanization  such  as  local  level,  space  stabilized, 
wander  azimuth,  free  azimuth,  strapdown,  etc. 

The  choice  of  navigation  error  variables  for  analysis 
has  often  followed  from  the  coordinate  system  used  in  the 
navigation  equations  or  that  implied  by  the  physical  I.N.S. 
platform  orientation.  Regardless  of  the  differences  in  the 
sets  of  equations  the -fact  is  that  I.N.S.  error  propagation 
is  to  a  large  extent,  completely  independent  of  system 
mechanization.  Britting  [Ref.  2]  has  shown  that  the  basic 
error  differential  equations  for  any  I.N.S.  may  be  written 
in  standard  coordinates,  regardless  of  the  physical 
nechani zation  or  internal  navigation  variables.  Further¬ 
more  ,  the  unforced  (homogeneous''  portion  of  these 
differential  equations  is,  under  certain  very  broad 
assumptions ,  identical  for  any  arbitrarily  configured 
terrestrial  I.N.S. 

The  error  equations  presented  in  this  chapter  follow  the 
philosophy  of  [Ref.  2]  including  the  choice  of  north-slaved 
coordinates  for  the  error  variables  and  the  identification 
of  the  unforced  (homogeneous)  portion  of  the  differential 
equations  that  is  independent  of  system  mechanization.  Two 
of  the  differences  that  are  noticeable  are: 

1)  The  error  equations  are  written  as  a  system  of  nine 
first  order  equations  (and  further  on  reduced  to  seven) 
rather  than  three  second  order  plus  three  first  order 
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equations.  The  first-order  form  is  the  state  space  repre¬ 
sentation  of  the  error  equations  used  in  modern  estimation 
theory . 

2)  The  form  of  altitude  compensation  assumed  in 
[Ref.  2]  is  not  found  in  most  inertial  navigators.  Gravity 
is  assumed  computed  as  a  function  of  the  inertial  system 
indicated  position. 

3.  GENERAL  I.N.S.  ASSUMPTIONS 

The  assumptions  pertaining  to  the  general  error 
differential  equations  are  broad  enough  to  encompass  all  of 
the  important  I.N.S.  configurations  [Ref.  2], 

1)  Three  accelerometers  are  available  to  measure  the 
specific  force  vector.  The  equations  for  a  two-accelero¬ 
meter  local  level  system  are  the  same  provided  the  inertial- 
altitude  and  vertical-velocity  equations  are  deleted. 

2)  The  accelerometer s  are  mounted  on  a  platform  whose 
angular  velocity  is  either  controlled  (as  with  gyro- 
stabilized  gimbaled  platforms)  or  is  measured  (as  with  the 
strapdown  systems)  . 

3)  The  system's  indicated  velocity  vector  and  three- 
dimensional  indicated  position  are  obtained  by  integration 
of  the  gravity  field  compensated  specific  force  measure¬ 
ments,  using  a  correct  set  of  differential  equations. 

4)  A  model  of  the  earth's  gravity  field  is  used  to 
calculate  the  gravity  vector  as  a  function  of  the  system- 
indie  a  t  *  d  po  s  i !:  •'  ;;  rt  . 
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5)  External  altitude  information  and  other  navigation 
measurements  may  be  used  to  update  the  inertial-navigation- 
system  indicated  position,  velocity,  and  attitude. 

6)  A  computer  is  availabe  to  process  the  navigation 
information  and  the  computation  errors  are  eitner  negligible 
or  may  be  treated  as  equivalent  instrument  uncertainties. 

?)  Both  the  mechanical  coordinate  frame  (the  frame 
tracked  by  the  platform)  and  the  computation  frame  (the 
frame  to  which  the  specific  force  measurements  are 
transformed  for  velocity  and  position  integration)  are 
arbitrar y . 

C.  LOCAL-LEVEL  TERRESTRIAL  NAVIGATOR 

Many  inertial  navigators  do  use  a  local-level  coordinate 
system  for  the  velocity  and  position  integration.  The 
local-level  terrestrial  navigation  system  physically  instru¬ 
ments  the  local  geographic  coordinate  frame.  The  platform 
axes  are  commanded  into  alignment  with  the  local  north-east- 
down  coordinate  system. 

The  local-level  terrestrial  system  is  undoubtedly  the 
most  successfuly  of  all  iner tial-n&vigation-system  configur¬ 
ations.  The  class  of  local-level  systems  today  constitutes 
the  majority  of  operational  inertial  navigations  systems. 
Since  this  system  is  described  in  detail  in  [Ref.  2]  here 
are  listed  some  of  its  advantages: 


1)  The  computation  of  gravity  components  is  greatly 
simplified.  In  fact  some  navigators  use  zero  for  both 
horizontal  components  of  gravity. 

2)  Some  inertial  navigators  have  no  vertical  accelero¬ 
meter  (this  is  the  case  in  the  present  study)  and  do  not 
mechanize  a  vertical  channel.  The  horizontal  velocity  and 
position  equations  of  a  local-level  set  are  appropriate  for 
such  a  nav igator  . 

3)  The  well-known  altitude  and  vertical  velocity 
instability  of  a  pure  inertial  navigator  must  be  stabilized 
by  means  of  an  external  altitude  reference.  But  a  local- 
level  set  of  variables  includes  altitude  and  vertical 
velocity  explicitly.  The  altitude  stabilization  equations 
therefore  can  be  simplified. 

4)  The  calculations  required  oo  provide  navigation 
outputs  and  displays  in  geographic  coordinates  are 
simplified . 

D.  THE  TWO  ACCELEROMETER  LOCAL-LEVEL  SYSTEM 

Many  inertial  navigators  have  only  two  accelerometers 
after  the  vertical  accelerometer  has  been  eliminated.  The 
system  is  composed  of  a  three  axis  inertial  platform,  two 
accelerometers  which  are  nominally  orthogonally  mounted  in 
the  instrumented  east  and  north  directions  and  a  computer 
which  performs  the  necessary  navigational  computations 
[Ref.  23. 
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The  north  and  east  gyros  are  respectively  connected  with 
the  instrumented  north  and  east  accelerometers  at  the  signal 
level  and  the  gyros  are  torqued  at  a  rate  proportional  to 
the  vehicle's  longitude  and  latitude  rates  so  that  the 
platform  can  maintain  its  axes  aligned  with  geographic  axes 
since  the  vehicle  carrying  the  navigation  system  is  assumed 
to  move  freely  over  and  above  the  earth.  The  accelerometer 
outputs  provide  these  required  torquing  signals  which  must 
be  so  compensated  that  gyro  command  can  be  obtained  as  a 
function  of  only  the  north  and  east  velocity  rates. 

Such  a  two-accelerometer  local-level  I.N.S.  has  seven 
state  variables:  two  of  position,  two  of  velocity,  and 
three  of  platform  alignment. 

1  •  Error  Model  aquations 

The  general  model  of  local-level  inertial  navigation 
systems  is  given  by  the  following  matrix  equations: 

Ax  =  ^  (20) 


where 

A  =  system  characteristic  matrix,  same  for  all  I.N.S. 
conf ig  urations 

g  =  forcing  vector  of  inertial  system  errors 

T 

or  s  £ 3 i »  ^2’  q3 '  ^4 •  ^5 »  ^6'^  (21) 
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x  s  error  state  vector  of  attitude,  position,  and 
velocity  errors 

T 

orj(  =  C£|^i6jr>6^,  S  L ,  6 1  ,  $  h  ]  (22) 

It  is  quite  important  to  emphasize  that  a  computer  simula¬ 
tion  program  developed  in  accordance  with  the  above  equation 
(20)  is  valid  for  all  possible  I.N.S.  configurations,  space 
stabilized,  strapdown  ,  wander-azimuth  and  not  only  for  the 
local-level  one.  Both  the  coordinate  frame  mechanized  by 
the  inertial  instruments  and  tne  computation  frame  are 
completely  arbitrary.  It  is  only  the  forcing  function,  q, 
which  depends  on  the  system  configuration  through  the  angu¬ 
lar  velocity  and  orientation  of  the  inertial  instruments. 

In  order  to  rewrite  this  equation  as  a  first  order 
vector-di f ferential  equation  tne  error  state  of  the  I.N.S. 
is  defined  as  [Ref.  2]: 

•t* 

e  (t)  =  [e^,  eg-,  c  q  ,  SL,  51  ,  SL,  51,  5h ,  5h]  (23) 

and  for  the  case  of  this  study  for  two-accelerometer  local- 
level  I.N.S.  system  this  reduces  to: 

e^(t)  =  [e^,  eg,  Eq,  SL,  Sl,  SL,  Si]  (24) 

where  the  seven  basic  I.N.S.  errors  are: 

eN'!  eE'  *D  =  north,  east,  down  platform  tilt  errors 

sL,  6 1  s  latitude,  longitude  position  errors 
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F6 2  =  *fD/r 


F63  =  fE/r 


^64  =  -1(1  +  2u>ie)cos  2L 


Fgy  =  -A  sin  2L 


=  fp/r  cos  L 


F73  =  “FN/r  GOS  L 


F Y4  —  1  u3n  L  +  2 A  L 


F?^  =  2x  tan  L 


F 77  =  2  L  tan  L 

i  ( 


£  s  7  x  1  error  vector  matrix 


wher« 


r  •  •  T 

.£  3  L  €  c  61,  5L,  6 1 J 

Ge  s  7  x  5  forcing  matrix,  with  non-zero  elements 


G11  =  1 


G22  =  1 


G33  =  1 
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5L,  ai  =  latitude,  longitude  rate  errors 

The  above  statements  allow  equation  (20)  to  be  written  as: 

e  =  Fg_e  +  Gcq  (25) 

where : 


12 


13 


14 


F 


34 


F 


37 


F 
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7x7  error  matrix  with  non-zero  elements 
- \  sin  L 

L 

-x  sin  L 

cos  L 

\  sin  L 

A  cos  L 

-1 

-L 

-x  cos  L 

-X  cos  L 

-  sin  L 

1 
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G64  =  1/r 


=  1/r  cos  t 


and  £=5x1  forcing  vector  matrix 
where 

a  =  q2»  q3*  aa,  q5lT  (23) 

and  neglecting  both  gyro  and  accelerometer  non-orthogonality 
errors  the  forcing  functions  are  comprised  of  10  I.N.S. 
component  errors  as  below: 


z 

3 

r _ 

"  "n  0  0  ' 

—  I  /■  1 1  N  |  r*  A 

!  P  j  ^u)cjh;  ;  +  "p 

0  rn  0 

1 

I  i(u)*d| 

o 

o 

~S 

j  L  - 

L  J 

( u)  f , 

i 

(u)f. 


where  for  our  study  of  two  accelerometer  local-level  system 
the  platform-to-navigation  transformation  matrix  Cpn  and  the 
I.N.S.  platform  torquing  rate  u»i p p  are  as  below: 


C  "  =  I 


x  cos  L 
-L 

-x  sin  L 
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TABLE  I 


I.N.S. 

SYSTEM 

MODEL 

MATRICES 

F  AND 

e 

G 

e 

0 

F 

r  12 

f13 

F 1  4 

0 

0 

F17 

F21 

0 

F23 

0 

0 

F26 

0 

:  3 1 

F32 

0 

?34 

0 

0 

u* 

1  37 

0 

0 

0 

0 

0 

F46 

0 

0 

0 

0 

0 

0 

0 

F57 

0 

F62 

F63 

c 

*  64 

0 

0 

F67 

c 

*71 

0 

f 

*  73 

‘74 

c 

k* 

‘  7o 

C* 

4  7  1 

i  ■ 

G11 

0 

0 

0 

0 

0 

G22 

0 

0 

0 

0 

0 

G33 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

°64  0 

0 

0 

0 

0 

G75 

(26) 


(27) 
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Finally  the  forcing  vector  matrix  is  modeled  for  our  case 


as  : 


r  i 

-  *“ 

rr 

(u)aN  +  rNX  cos  L 

*2 

(  u  ) 

<53 

- 

(u)u»D  -  rQx  sin  L 

(u)fN  -  5g 

q5 

(U)  f  r-  +  n , 

L  a 

2 .  Error  Equations  Solution 

The  solution  of  the  differential  equations 
represented  by  equation  (25)  gives  the  error  response  for 
the  two-accelerometer  local-level  navigator  for  arbitrary 
vehicle  motion  within  the  constraints  implied  by  a  "first- 
order"  analysis.  Since  the  coefficients  of  the  differential 
equations  are  time  varying  the  analytic  solution  of  equation 
(25)  would  be  quite  tedious  and  require  the  user  of  a 
computer  program  to  generate  flight  profiles.  In  our  study 
specific  cases  are  examined  so  that  the  coefficients  of  the 
differential  equations  can  easily  be  calculated. 

The  specific  cases  which  are  examined  are  the 
following  three: 

••  ••  •• 

•  •  •  • 

1)  Stationary  case,  where  x  =  and  L=L=l=l=h=h  =  0 

2)  Easterly  flight  at  600  "ft/sec"  or  355-5  "knots," 
where  x  =  1.557  x  u>ie 


3)  Westerly  flight  at  600  ,,ft/sec"  or  355.5  "knots," 

where  x  =  0.442  x  «. 

le 

Writing  the  equation  (25)  in  terras  of  the  error 
states  explicitly  we  have  the  following  system  of  differen¬ 
tial  equations  which  must  be  solved  simultaneously: 

x(1)  =  -x(sinL)x(2)  -  x(sinL)x(4)  +  (cosL)x(7)  +  q  ■] 

x ( 2  )  =  x ( sin  L ) x ( 1 )  +  x(eosL)x(3)  -  x(6)  +  q5 

x(3)  =  -x(cosL)x(2)  -  X(cosL)x(4)  -  (sinL)x(7)  + 

x(4)  =  x ( 6 )  (32) 

x (5 )  =  x (7 ) 

x  (  6  )  =  =  x  ( 2)  -  x  ( s i n 2 L  )  :<  (  7  )  +  q4 

x(7)  =  -  r (  coslT  x(1)  +  2x ( tan L )  x (  6  )  +  q5 

The  values  of  the  parameters  used  for  the  computer 
simulation  are  given  in  Table  II. 

a.  Constant  Gyro  Drift  Errors 

Letting  constant  gyro  drift  be  the  sole  error 
source  in  the  I.N.S.  system  where  the  constant  gyro  drift 
rates  (u)w^,  (u)uj.,  (u)a)^  are  associated  with  the  north, 
east,  and  azimutn  gyros  respectively,  computer  simulation 
verfied  the  following. 
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TABLE  II 

COMPUTER  SIMULATION  PARAMETERS 


Geographic  Latitude,  L 

45°L 

Terrestrial  Longitude,  1 

0°  (Greenwich) 

Earth  Rate,  <uie 

Constant  Gyro  Drift 

0.2618  Rad/hour 

•  Nortn,  (u)^ 

1  meru  ( 1 ) 

•  East,  ( u ) y - 

1  meru 

•  Azimuth,  (u)uiq 

1  meru 

Constant  Accelerometer  3ias 

•  North,  (u)fN 

200  ug  (2) 

•  East,  (u)f£ 

200  ug 

Initial  Platform  Tilt 

•  North,  6^(0) 

0.14  mrad 

•  East,  e- ( 0  ) 

0.14  mrad 

•  Azimuth,  e^(0) 

0.14  mra  i 

Initial  Latitude  Error,  5L(0) 

0.17  mrad  (3) 

Initial  Longitude  Error,  61(0) 

0.17  mrad 

Initial  Latitude  Rate  Error,  6L(0) 

0.34  mrad  (4) 

Initial  Longitude  Rata  Error,  51(0) 

0 . 34  mrad 

Stationary  Flight 

)  =  "ie 

Constant  Easterly  Flight  at  600  ft/sec 

x  =  1  •  557xai. 

1  e 

Constant  Westerly  Flight  at  600  ft/sec 

X  =  0.442xaie 

Gravitational  Acceleration  Constant 

g  =  32. 2  ft/sec 

Earth  Radius 

r  =  6,37  8-3  Km 

1.  The  RMS  gyRO  drift  of  1  meru  corresponds  to  0.015°/'ir  or 
0.2618x10  3  rad/hr. 

2.  The  RMS  accele 'ometec  bias  of  200ug  corresponds  to 

0. 4356x  1 0”-5  r  q  j/  ( hr  )  . 

3.  The  RMS  position  error  of  0.17  mrad  corresponds  to  1085m 
or  0.586  arc  min. 

4.  The  RMS  velocity  error  of  0.34  mrad  corresponds  to 

2  ft/ sec. 
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For  the  stationary  case  we  observe  in  Figures  4 
through  10  that  for  the  north  and  east  level  errors,  and 
eg,  the  Foucault  modulation  is  an  effect  of  first-order  in 
contrast  witn  the  latitude,  longitude  and  azimuth  errors  6 L, 
si,  and  eD  respectively  where  the  Foucault  modulation  has 
only  a  second-order  effect  [Ref.  1].  These  computer 
solution  results  suggest  that  it  would  be  convenient,  for 
design  purposes,  to  neglect  tne  Foucault  modulation  since 
tne  equations  we  obtain  then  are  easily  solved  and  give 
solutions  with  approximately  the  same  amplitude  information 
for  latitude  and  longitude  which  are  of  primary  importance 
for  navigational  purposes  while  the  relatively  poor 
information  of  the  level  errors  is  of  secondary  importance. 

We  further  note  from  these  computer  graph 
solutions  that  tne  effect  of  tne  Foucault  cerms  in  tne  error 
equations  system,  is  to  modulate  the  Schuler  oscillations  at 
a  frequency  given  by  the  local  vertical  projection  of  earth 
rate,  namely  <d,  x  sin  L,  which  corresponds  to  a  period  of 

Q 

33-9  hours  for  the  selected  latitude  L  =  45  • 

For  the  case  of  constant  easterly  flight  at  600 
ft/3ec  the  results  are  given  in  Figures  11  to  17.  Compari¬ 
son  with  the  curves  for  the  stationary  case  indicates  that 
the  lowest  modulation  frequency  has  increased  from  x  =  <"^e 
to  x  s  1.557  x  u)ie  and  the  space  rate  period  is  10.9  hours 
while  the  Foucault  modulation  now  occurs  with  a  period  of 
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about  21.8  hours  instead  of  the  33*9  hours  period  for  the 
stationary  case. 

Another  important  feature  revealed  by  the  above 
comparison  is  that  the  azimuth  and  latitude  errors  are 
reduced  from  the  corresponding  for  the  stationary  case  by  a 
factor  of  1.557  which  represents  the  ratio  for  this 

case  . 

For  the  responses  to  the  north  and  azimutn  zero 
drift,  (u)uu  and  (u)u>r,  the  vehicle  motion  appears  to  have 
little  effect  on  the  error  growth  in  the  cases  that  exhibit 
a  longitude  error  which  grows  with  time. 

The  level  errors  in  response  to  level  gyro  drift 
are  seen  to  remain  unchanged  while  the  level  error  response 
to  azimuth  gyro  drift,  (u)«-y  is  seen  to  emerge  from  the 
computer  noise  having  a  peak  value  of  2.3  rad/meru  (Figure 
13) . 

The  longitude  error  in  response  to  azimuth  gyro 
drift  which  was  bounded  for  the  stationary  case  is  now 
reduced  by  the  factor  x/(uie  or  by  1.55  while  the  latitude 
and  longitude  rate  error  magnitudes  are  unaffected  by  the 
vehicle  motion . 

Finally  for  the  westerly  flight  case  it  is 
verified  that  the  level  errors  remain  unchanged  without  the 
effect  of  Foucault  modulation,  but  the  latitude,  longitude, 
and  azimuth  errors  grow  approximately  in  proportion  to  the 
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time-drift  rate  prodct.  The  computer  solution  graphs  for 
the  westerly  flight  case  are  shown  in  Figures  19  through  24. 

Similar  results  are  found  for  the  cases  of  east 
and  azimuth  gyro  drift  but  they  are  not  included  here  due  to 
large  amount  of  graphs. 

b.  Accelerometer  Bias  Errors 

Considering  the  accelerometer  bias  as  the  sole 
error  source  computer  simulation  of  equations  ^25)  shows  the 
following  result  of  the  effects  of  the  north  and  east 
accelerometer  bias,  (u)fN  and  (u)f£  respectively,  on  the 
navigation  and  level  errors. 

For  the  stationary  case  the  results  are  shown  in 
Figures  25  through  3 1 •  We  note  that  the  Schuler  mode  is 
predominant  since  the  accelerometer  bias  directly  excites 
the  relatively  high  gain  level  loops  and  that  the  Schuler 
oscillations  are  modulated  at  the  Foucault  mode  frequency  of 
33.9  hours  per  cycle.  The  maximum  values  for  the  navigation 
errors  proved  to  be  in  the  range  of  20  x  10  rad/200  ug  for 

the  latitude  error  and  1.4  x  10“^  rad/200  ug  for  the 
longitude  error. 

We  notice  as  for  the  constant  gyro  drift  case 
that  we  can  neglect  the  effects  of  Foucault  modulations  as 
first-order  ones  and  proceed  in  the  solution  of  the 
resulting  equations  more  easily  obtaining  almost  the  same 
approximate  amplitude  information  of  the  navigation  and 
level  errors. 
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For  the  case  of  constant  easterly  flight  at  600 
ft/sec  the  computer  simulation  graphs  are  shown  in  Figures 
32  through  38.  We  observe  here  that  again  the  Foucault 
modulating  frequency  has  increased  by  a  factor  of  1.55  which 
corresponds  to  the  ratio  x/uiig.  Nevertheless  we  see  that 
the  error  sensitivities  remain  unchanged  with  the  previously 
explained  stationary  case. 

Figures  39-45  show  computer  solutions  of  the 
navigation  and  level  errors  for  the  case  of  westerly  flight 
at  600  ft/sec  . 

We  see  now  that  the  Foucault  modulating 
frequency  nas  decreased  oy  a  factor  0.442  corresponding  to 
the  ratio  \/  in  this  case  and  once  again  we  can  proceed 
to  the  solution  of  the  equations  without  considering  the 
Foucault  terms,  especially  for  design  purposes.  An  easy 
extension  of  the  above  observations  is  that  for  the  limiting 
case  when  the  terrestrial  longtitude  rate,  i,  in  a  western 
direction  is  equal  to  the  earth  rate,  u>ig,  then  the  Foucault 
modulation  disappears  completely  leaving  a  pure  Schuler 
osc illation . 

c.  Initial  Condition  Errors 

The  results  on  the  navigation  and  level  errors 
due  to  effects  of  the  initial  conditions  are  now  presented 
accompanied  by  only  the  most  important  graphs  of  the 
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computer  simulation  since  the  total  amount  of  graphs  is  too 
large  to  be  included  in  this  study. 

For  the  stationary  case  selected  graphs 
representing  the  level  errors  for  an  initial  north,  east  and 
azimuth  level  error  of  0.14  mrad  or  0.438  arc-min  are  shown 
in  Figures  46-47,  43-49,  and  50-51  respectively. 

Figures  52  through  57  present  the  level  and 
navigation  errors  for  an  initial  latitude  rate  error  of  2 
ft/sec  corresponding  to  0.34  mrad/ hour  while  figures  53 
through  63  show  the  associated  errors  with  an  initial 
longitude  rate  error  of  the  same  amount. 

There  is  no  need  to  include  any  graphs  for  the 
resulting  errors  due  to  initial  longitude  error  since  by 
inspection  of  the  error  differential  equations  we  can  see 
that  longitude  is  uncoupled  from  the  other  computation  loops 
so  that  an  initial  longitude  error  holds  constant  and  no 
other  error  becomes  non-zero. 

We  discussed  up  to  now  the  errors  of  a  pure 
I.N.S.  system.  In  the  next  chapter  we  proceed  with  the 
consideration  of  the  combined  I. N.S./3. P. S.  system  and  the 
results  we  achieved  after  computer  simulation. 


IV.  I.N.S./G.P.S.  SYSTEM  MODEL  AND  EQUATIONS  SOLUTION 


In  order  to  apply  the  Kalman  filter  equations  discussed 
in  Chapter  II  a  reference  system  model  which  is  a  good 
approximation  to  the  real  world  dynamics  is  needed. 

In  this  chapter  we  outline  the  reference  I. N. S./G . P. S. 
system  equations  selected  for  this  study.  First  we  are 
defining  the  error  states  incorporated  in  the  system  model 
along  with  their  assumed  initial  conditions  values.  Next  we 
discuss  the  modeling  of  the  I.N.S.  plant  error  states  and 
finally  we  present  the  equations  for  the  integrated 
I.N.S./G.P.S.  system  together  with  their  simulated  computer 
results . 


A.  SYSTEM  MODEL 

1 .  State  Variable  Definition  and  Initial  Conditions 
In  Table  III  we  present  a  listing  of  the  state 
variables  utilized  in  the  reference  system  model.  The 
initial  conditions  on  the  I.N.S.  error  states  are  highly 
arbitrary  and  the  selected  values  are  similar  to  those  used 
in  other  unclassified  studes  [Refs.  4,5]. 

For  the  initial  conditions  on  the  gyro  error  states 
and  the  accelerometer,  the  values  are  selected  for  a  typical 
inertial  navigation  system  of  one  nautical  mile  per  hour 


After  the  above  definition  of  the  initial  conditions 
we  have  by  the  same  time  specified  in  a  complete  way  the 
initial  covariance  matrix  P(0),  since  its  diagonal  elements 
are  the  squared  values  of  the  given  RMS  initial  conditions. 
The  remaining  off-diagonal  elements  of  the  initial  covar¬ 
iance  matrix  are  assumed  to  be  zero  initially. 

Furthermore  the  propagation  of  the  linear  variance 
equation  (13)  requires  an  additional  knowledge  of  the  two 
matrices  F  and  Q*  where: 

Q*  =  GQGT  (33) 

where 

G  is  the  forcing  input  matrix 

Q  is  the  input  noise  covariance  matrix 

and  the  F  matrix  is  the  same  as  in  Equation  (25)  and  which 
has  been  used  in  the  previous  chapter  for  the  inertial 
navigation  system  error  equations  solution  and  computer 
simulation . 

For  the  Q*  matrix  the  only  non-zero  elements  are  all 
diagonal  and  we  will  denote  these  from  now  on  as  where 
the  subscript  i  denotes  the  row  and  column  of  the  value. 

For  example,  indicates  that  this  is  the  value  which 
belongs  to  the  intersection  of  the  3rd  row  and  the  3 r*d 


column  in  the  Q*  matrix  and  corresponds  to  a  white  noise 
input  on  state  variable  number  3. 

These  non-zero  elements  in  the  reference  system  Q* 
matrix  are  five,  corresponding  to  3tate  numbers  1,  2,  3»  6 
and  7  according  to  the  notation  of  Table  III. 

2.  P_Lant  Error  States 

The  following  seven  states,  North,  East,  and  Azimuth 
level  errors,  X  and  Y  position  errors,  X  and  Y  velocity 
errors,  constitute  the  plant  error  states.  The  differential 
equations  of  these  states  describe  the  natural  unforced 
dynamic  response  of  the  errors  in  the  inertial  navigation 
system . 

Tnere  are  various  models  for  the  implementation  of 
these  error  states.  As  we  did  in  the  previous  chapter  we 
will  use  3gain  the  Pinson  error  model  descrioed  by  the 
matrix  F  given  in  equation  (25)  for  our  specific  case  of  the 
local-level  two- accelerometer  inertial  navigation  system 
conf ig uration  . 


B.  EQUATIONS  SOLUTION  AND  COMPUTER  SIMULATION 

Using  the  definitions  described  in  previous  pages  we  can 
write  the  following  equations  for  the  error  states: 

x(  t)  =  F  x(  t)  +  G  w(t) 

(34) 

z(  t)  s  H  x(t)  +  v(t) 
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TABLE  III 


I.N.S./G.P.S.  SYSTEM  STATE  VECTOR  DEFINITION 

Error 


State 

Symbol 

Definition  RMS  Initial 

Condition 

I.N.S.  PLANT  ERROR  STATES 

1 . 

e  N 

North  atttiude  error 

0.  14x10“1 2 3 

Rad 

2 . 

e  „ 

U 

East  attitude  error 

0.  14x 10“3 * 

Rad 

3- 

4. 

6  D 

5L 

Azimuth  attitude  error 

Y  position  error 

0.  14x 10~3 
0.  17x 10'3 

Rad 

Rad 

( 1 ) 

5. 

6. 
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6  L 

X  position  error 

Y  velocity  error 

0. 1 7x 1 0~3 
0. 34x 10'3 

Rad 

(2) 

Rad/hour 

7  • 

d 

X  velocity  error 

0. 34x 10~3 

Rad 

•'hour 

I.N.S.  ERROR  SOURCES 

3. 

(u)  JN 

North  Gyro  Drift 

1  meru 

(3) 

9. 

(  U  )ojg 

East  Gyro  Drift 

1  meru 

10. 

(u)uD 

Azimuth  Gyro  Drift 

1  meru 

1 1 . 

(u)fN 

North  Accelerometer  Bias 

200x10~6g 

(4) 

12. 

(u)fE 

East  Accelerometer  Bias 

£ 

200x 10  g 

1.  The  RMS  position  error  of  0.17  milliradians  correponds  to 
1085m  or  0.586  arc  min. 

2.  The  RMS  velocity  error  of  0.34  rnillirad/hour  corresponds 
to  2  ft/sec. 

3.  The  RMS  zvro  drift  of  1  meru  corresponds  to  0.015°/hr  or 

261.8x10  rad/hour.  , 

4.  The  RMS  acceleroraetec  bias  of  200x10”°g  corresponds  to 

0. 4356x  10--3  rad/(hr)  . 
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where 


x(t)  =7x1  error  state  vector 
=  [eN.  e  £»  ®d*  6L’ 

F  =7x7  Pinson  error  model  matrix 
as  described  in  equation  (25) 

G  =  7  x  5  input  forcing  matrix  as  descrioed 
in  equation  (25) 

£(t)  =2x1  vector  of  measured  states 

H  =2x7  measurement  matrix  where 

[  0  0  0  1  0  0  0 

H  =  f 

10  0  0  0  10  3 

J- 

denoting  that  we  have  available  measured  information  for  the 
X  and  Y  position  error. 

w(t)  =5x1  forcing  vector  assumed  to  be  white 
Gaussian  noise 


and 


v(t)  =2x1  vector  of  measurements  noise  assumed  to 
be  white  Gaussian. 

Using  the  feedback  configuration  of  the  Kalman  filter  we 
can  write  the  following  equations: 
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(35) 


x(t)  =  F  x(t)  +  G  w ( t )  +  K[z  -  H  x(t) ] 

and 

P(t)  =  F  P  +  P  FT  +  G  Q  GT  -  P  HT  R"1  H  P 
where 

£(t)  =  the  covariance  matrix 
T  - 1 

-  =  L  H  1  the  Kalman  filter  gains  matrix 

£  =  the  measurement  noise  covariance  matrix 

In  order  to  achieve  numerical  results  via  computer 
simulation  we  write  the  predicted  error  states  of  our  system 
in  explicit  form  as  below  with  the  help  of  Table  IV  in  which 
the  system  states  and  their  corresponding  symbols  are 
defined : 

x(1)  =  -i(sinL)x(2)  -  i(sinL)x(4)  +  (cosL)x(7)  +  AA  + 

♦  kUCx(8)  -  x  ( 4)  ]  +  K12[x(9)  -  x  (5)  ]  (35-a) 

x(2)  =  A(sinL)x(l)  +  x(cosL)x(3)  -  x(6)  +  BB  + 

+  K2iCx(8)  -  x ( 4 ) ]  +  K22[x(9)  -  x (5) 3  (35-o) 

a 

x(3)  =  -i(cosL)x(2)  -  i(cosL)x(4)  -  (sinL)x(7)  +  SS  + 

*  K3 T Cx ( 8)  -  x ( 4 ) ]  +  K32[x(9)  -  x (5) 3  (35-c) 

x(4)  s  8(6)  +  K41C8(8)  -  8(4)3  ♦  Ki42Cx(9)  -  8(5)3  (35-d) 


64 


TABLE  IV 


COMPUTER  SIMULATION  VARIABLES  AND  CONSTANTS 


x(  1) 

z 

e  N 

z 

North  attitude  error 

x  (2) 

- 

e  E 

z 

East  attitude  error 

X C  3) 

z 

gd 

z 

Azimuth  attitude  error 

X  ( 4 ) 

z 

«L 

= 

Y  position  error 

x(5) 

- 

U 

z 

X  position  error 

x  (6 ) 

z 

6  L 

z 

Y  velocity  error 

x  ( 7 ) 

= 

si 

z 

X  velocity  error 

x  ( 3 ) 

= 

SLg 

= 

G.P.S.  Y  position  error  measurement 

x  ( 9 ) 

z 

513 

z 

G.P.S.  X  position  error  measurement 

x(  10) 

z 

5Lt 

= 

True  Y  position  error 

x(  1  1) 

z 

slt 

z 

True  X  position  error 

x  (  12) 

= 

w 

z 

Input  white  Gaussian  noise 

x(  13) 

z 

V 

z 

Measurement  white  Gaussian  noise 

A 

z 

[(u)o»N]2 

z 

North  Gyro  Drift  Variance 

B 

= 

[(u)o>E]2 

z 

East  Gyro  Drift  Variance 

S 

z 

C(u)u.D]2 

= 

Azimuth  Gyro  Drift  Variance 

D 

z 

c (u) f n32 

z 

North  Accelerometer  Bias  Variance 

E 

z 

[(u)fE]2 

z 

East  Accelerometer  Bias  Variance 

F 

z 

R11 

z 

Measured  Y  position  error  variance 

G 

z 

R22 

= 

Measured  X  position  error  variance 

AA 

z 

x( 12)*A 

s 

White  noise  like  North  Gyro  Drift  strength 

BB 

z 

x(  12) *  B 

z 

White  noise  like  East  Gyro  Drift  strength 
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TABLE  IV  (CONTINUED) 


=  White  noise  like  Azimuth  Gyro  Drift 
strength 

a  White  noise  like  North  Accelerometer 
bias  strength 

=  White  noise  Like  East  Accelerometer 
bias  strength 

a  White  noise  like  Y-position  error  strength 
=  White  noise  like  X-position  error  strength 
=  constant 


x (5)  =  x (7 )  +  K51[x(8)  -  X ( 4 )  3  +  K52[x(9)  -  5(5)]  (35-e) 

x ( 6 )  =  kx(2)  -  x ( sin2L) x ( 7 )  +  DD  + 

+  K6l[x(8)  -  x  (4 )  ]  +  K62[x(9)  -  5(5)]  (35-f) 

x  (7)  =  -  *(D  +  2x(tanL)x(6)  +  EE  + 

COSL 

+  K?1  [  x  (  8)  -  5(4)3  +  K?2[x(9)  -  5(5)]  (35-g) 

Assuming  the  input  forcing  vector  as  white  Gaussian 
noise  whose  strength  is  related  to  the  value  of  the  variance 
of  each  input  error  source  (the  corresponding  one)  computer 
simulation  was  proceeded  in  the  following  way. 

First  with  the  help  of  the  RICATI  FILTER  computer 
program  available  at  the  NFS  *J.R.  Church  Computer  Center  we 
solved  the  cor  responding  for  our  study  Ricati  equation  of 
covariance  propagation  obtaining  the  Kalman  filter  gain 
matrix.  The  data  we  used  to  run  the  above  program  are 
provided  in  Tables  III  and  V.  A  listing  of  the  data 
formulation  for  the  RICATI  FILTER  program  is  given  in 
Appendix  D. 

The  calculated  from  the  above  program  values  of  the 
Kalman  filter  gains  for  a  processing  period  of  four  hours 
are  given  in  the  following  Table  VI. 

Additional  runs  of  the  above  program  have  been  contacted 
for  processing  periods  up  to  36  hours  and  it  has  been 
observed  that  the  Kalman  filter  gains  reach  a  steady  state 


67 


TABLE  V 


NUMERICAL  VALUES  FOR  RICATI  PROGRAM1 
F  Matrix  (7x7) 


12 

= 

-0.1851 

F„  = 

-0.  1851 

F 17  = 

0.707 

21 

s 

0.1851 

F23  ■ 

0.1851 

F26  = 

-1.0 

32 

= 

-0.  1851 

F34  = 

-0.  1851 

F37  = 

-0.707 

46 

= 

1.0 

57 

= 

1 . 0 

62 

s 

19-92 

F67  = 

-0. 2618 

71 

= 

-28. 175 

F76  - 

0.5235 

All  other  elements  are  zero. 


G1  Matrix  (5x7) 

!j11  =  a2  2  =  G33  =  1,0 
G46  =  0.0915 

G57  =  0. 1294 

All  other  elements  are  zero . 


H  Matrix  (2  x  7) 


H 1 4  =  H25 


1.0 


All  other  elements  are  zero. 


For  the  used  values  of  elements  in  the  matrices  the 
results  will  be  given  having  units  the  appropriate  for  each 
case  combination  of  radians  and  hours. 
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TABLE  V  (CONTINUED) 


Q  Matrix  (5  x  5) 

Qn  =  Q22  =  Q33  =  [0.0002618  rad/hr]2  =  0. 7x10'?[rad/hr]2 
Q44  S  Q55  =  [0.0004355  rad/(hr)2]2  =  0. 19x10"6[rad/(hr)2]2 

All  other  elements  are  zero. 


R  Matrix  (2  x  2 ) 


R11  =  R22 
R 1 2  =  R21 


/*  p 

0.52x10'°  [rad]“ 

0.0 


?(0)  Matrix  (7  x  7) 


Pn(3)  =  ?p2  (  3  )  =  ?  3  3  ( 0  >  =  [0.00014  rad]2  =  0  .  1  xl  0‘  7  [  rad  ]  2 
P44(0)  =  ?55(0)  =  [0.00017  rad/nr]2  =  0 . 3 xl 0" 7 [ r ad/hr ] 2 
P66(0)  =  P77 < 0 )  =  [0.00034  rad/(hr)2]2  = 

0.115  x  10’6  [rad/ (hr)2]2 

All  other  off-diagonal  elements  are  zero. 
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table  VI 


KALMAN  FILTER  CASE  FOR  A  4  HOUR  PROCESS 


K11 

= 

0.034756672 

K12 

=  1.17201123 

*21 

= 

-1 . 36921265 

<22 

=  -0.0401669175 

K3i 

= 

1.56555360 

<32 

=  -2.86193609 

K41 

= 

3. 4402061 1 

<42 

=  0.090728150 

K51 

s 

0.0907281520 

<52 

=  4.31037639 

<61 

= 

5.92162434 

• 

<52 

=  -0.339257459 

K71 

= 

1 . 04245351 

<72 

=  9.29378811 
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condition  for  which  their  values  are  not  very  much  different 
from  those  achieved  for  a  4  hour  processing  period. 

So  in  the  following  calculations  and  computer  simulation 
we  have  used  the  values  of  the  Kalman  filter  gains  resulted 
for  a  4  hour  simulation  period. 

Having  available  the  values  of  Kalman  filter  gain  matrix 
elements  which  are  used  to  multiply  the  residuals  in  the 
appropriate  equations  in  order  to  achieve  the  predicted 
error  states  of  the  integrated  I. M.S./G. P.S.  system,  the 
appropriate  program  has  been  formulated  in  order  to  solve 
the  error  differential  equations  (32)  described  above,  with 
the  use  of  the  available  routine  INTEG2S  slightly  modified 
for  accurate  evaluation  and  plot  of  the  error  state 
variables . 

The  simulation  results  for  the  I.N.S./G. P.S.  system 
operation  for  a  period  of  four  (4)  hours  are  presented  below 
in  Figures  64  through  93-  We  can  easily  observe  that  all 
the  state  error  variables  of  the  combined  I.N.S./G.P.S. 
system  are  damped  out  and  the  resulting  value  of  the  errors 
after  a  period  of  one  (1)  hour  is  small  enough  so  that  the 
operation  of  our  system  model  can  be  characterized  as 
successful . 

Specifically  in  Figure  64  we  observe  that  the  north 
level  error  of  our  system  is  dropped  down  to  0.025  milli- 
radians  after  one  (1)  hour  even  if  the  starting  initial 
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condition  value  was  0.14  railliradians .  Furthermore  at  the 
end  of  a  4  hour  period  the  error  has  been  diminished  to  the 
value  of  0.005  milliradians  which  is  subject  to  error 
reduction  by  a  factor  of  28. 

In  Figures  65  and  66  the  East  and  Azimuth  attitude 
errors  are  presented  respectively  where  similar  as  with  the 
north  attitude  error  observations  occur. 

The  Y  position  error  behavior  is  presented  in  Figure  67- 
There  we  can  see  that  even  if  we  started  from  an  initial 
condition  error  of  0.17  milliradians  (or  3256  ft)  after  one 
hour  processing  the  error  has  been  diminished  to  only  0.028 
milliradians  (or  536  ft)  and  furthermore  after  a  four  hour 
period  this  error  drops  down  to  0.01  milliradians  (or 
191.5  ft)  . 

The  X  position  error  damping  out  seems  to  be  more 
attractive  since  from  Figure  68  we  can  see  that  after  one 
hour  the  error  drops  down  to  0.020  railliradians  (or  383  ft) 
and  at  the  end  of  a  four  hour  period  the  error  is  diminished 
to  0.0002  railliradians  (or  3.83  ft)  which  is  very  small  con¬ 
sidering  also  that  we  started  with  an  initial  condition 
value  of  the  X  position  error  of  0.17  milliradians  (or 
3256  ft). 

Figure  69  presents  the  propagation  of  Y  velocity  error. 

It  is  easily  observed  that  this  error  drops  down  to  the  very 
small  value  of  0.040  tnilliradians/hour  (or  68  x  1 0-^  ft/sec) 
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after  one  hour  and  to  the  negligible  error  of  0.002 
milliradians/hour  (or  34  x  10"^  ft/sec)  which  provides  the 
advantage  of  very  accurate  evaluation  and  tracking  of  the  Y 
velocity  state  variable. 

Similar  with  the  above  considerations  and  even  better 
results  occur  for  the  case  of  the  X  velocity  error  of  the 
integrated  I , N. S. /G. P.S.  system.  We  see  from  Figure  70  that 
this  error  starting  from  an  initial  condition  value  of  0.34 
milliradians/hour  (or  2  ft/sec)  drops  down  to  0.1  milli¬ 
radians/hour  (or  17  x  10  ^  ft/sec)  after  one  hour  operation 
and  furthermore  down  to  0.002  milliradians/hour  (or  34  x 

_5 

10  ft/sec)  after  a  period  of  four  hours  which  again 
denotes  a  very  accurate  tracking  of  the  X  velocity  error 
state  variable. 

In  Figures  71  and  72  the  normalized  inserted  and 
measurement  noise  of  the  combined  I.N.S./G.P.S.  systems  are 
presented  respectively. 

Thinking  of  the  operation  of  our  system  in  the  long 
term,  results  of  the  computer  simulation  are  presented  in 
Figures  73  through  79.  Using  the  same  input  data  for  our 
I.N.S./G.P.S*  system  model  and  running  the  program  for  a 
36  hour  process  we  see  that  the  behavior  of  the  feedback 
configration  of  the  Kalman  filter  in  our  system  continues  to 
be  attractive  throughout  the  long  term  period  of  interest 
without  diverging  at  any  moment. 
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In  addition  to  the  above  considerations ,  in  order  to 
make  our  system  more  realistic  and  compatible  to  the  real 
world’s  conditions  we  put  some  noise  in  the  two  error  state 
equations  of  X  and  Y  position  which  did  not  include  any 
noise  from  our  theoretical  design  of  the  system.  So  we 
replace  the  two  equations  (35-d)  and  (35-e)  in  our  system  of 
equations  with  the  following  two  equations: 

X  (4)  =  8(6)  +  AA  K41[x(8)  -  x  (4)  ]  *  K42[x?9)  -  8(5)]  ( 35-d  ’  ) 
and 

x(5)  =  x ( 7 )  +  AA  +  K51[x(8)  -  8(4)]  +  K52Cx(9)  -  8(5)]  (35-e’) 

Assigning  to  the  strength  of  this  intentionally  inserted 
noise  a  value  similar  to  that  of  the  strength  of  the  gyro 
drift  (that  is  a  value  of  0.0685  x  10“°  [rad]-)  we  ran  the 
same  program  and  we  achieved  results  proving  that  the 
combined  I.N.S./G.P.S.  system  reacted  in  a  way  exactly  the 
same  as  it  had  reacted  without  the  inserted  noise  in  the  X 
and  Y  position  error  equations.  So  we  make  the  conclusion 
that  the  intentionally  inserted  noise  did  not  affect  the 
operation  of  our  system  model  neither  from  the  accuracy 
point  of  view  nor  from  the  time  point  of  view. 

The  above  considerations  and  results  can  be  seen  in 
Figures  80  through  86  for  the  four  (4)  hours  short  term 
process  and  in  Figures  87  through  93  for  the  36  hours  long 
term  operation. 
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In  Table  VII  on  the  next  page  we  summarized  the  state 

errors  of  the  combined  I. N. S./G. P. S.  system  after  a  period 

of  two  and  four  hours  operation.  In  the  same  table  we 

included  the  starting  initial  condition  for  each  error  state 

in  order  to  make  our  comparisons  easier  and  handy.  The 

results  included  in  Table  VII  are  those  achieved  from  the 

computer  simulation  without  any  noise  corrupting  the  two 

error  states  of  Y  and  X  position.  But  since  the  addition  of 

noise  with  strength  similar  to  that  of  the  gyro  drift 
6  2 

(0.0685  x  10  [rad]  )  di  not  affect  the  system  model 
operation  as  mentioned  before,  the  same  Table  VII  represents 
also  the  summary  of  state  errors  for  the  real  world's  system 
model  of  the  I. N. S. /G. P. S.  system. 

Up  to  now  we  considered  our  system  to  be  corrupted  by 
white  Gaussian  input  noise.  Since  in  the  real  world  in  many 
cases  the  presence  of  colored  noise  is  apparent  we  must 
consider  the  operation  of  our  I. N.S./G. P. S.  system  under  the 
presence  of  such  noise  and  compare  the  results  with  those 
achieved  when  the  system  was  driven  by  white  noise. 

In  the  following  section  a  realistic  modeling  of  the 
I.N.S.  component  errors  is  discussed  and  the  results  of  the 
computer  simulation  are  presented  together  with  the  compari¬ 
son  conclusions  of  the  system's  operation  under  colored 
noise  versus  white  noise  corruption. 
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TABLE  VII 


SUMMARY  OF  STATE  ERRORS  FOR  THE  I .N .S ./G. P.S.  SYSTEM  MODEL 


State 

Error 

Initial 

Condition 

Error  in  2  hours 

Error  in  4  hours 

e  N 

0.14  mrad 

0.025  mrad 

0.005  mrad 

eE 

0.14  rarad 

0.021  mrad 

0.001  rarad 

eD 

0. 14  mrad 

0. 04  mrad 

0.02  mrad 

5  L 

0. 17  mrad 
(3256  ft) 

0.028  mrad 
(536  ft) 

0.01  mrad 
(191*5  ft) 

6  I 

0.17  mrad 
(3256  ft) 

0.020  mrad 
(383  fft) 

0. 0002  mrad 
(3-3  ft) 

0. 34  mrad/ hr 
(2  ft/ sec ) 

0 . 040  mrad/ hr 
(0.0068  ft/sec) 

0. 002  mrad/hr 
(0. 00034  ft/sec) 

0.  34  mrad/ hr 
(2  ft/sec) 

0 . 1  rarad /hr 
(0.017  ft/sec) 

0.002  mrad/hr 
(0. 00034  ft/sec) 
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c.  I.N.S.  COMPONENT  ERROR  MODELS 

In  Chapter  III  equation  (3D  indicates  that  the  I.N.S. 
component  errors  consist  of  three  gyro  drift  uncertainties, 
two  accelerometer  measurement  uncertainties,  two  gyro 
torquer  scale  factor  errors  and  two  geodetic  uncertainties. 
Realistic  modeling  of  the  two  major  error  components,  the 
gyro  drift  and  the  accelerometer  measurement  is  described 
below . 

1 .  Gyro  Drift  Uncertainties 

The  three  gyro  drift  uncertainties,  (u)u>N,  (u)u>E, 
(u)u>D  are  each  modeled  as  an  exponentially-correlated 
(colored)  noise  plus  an  additive  random  (white)  noise: 

(u)'j|i  :  5.  +  w  ;  i  =  N ,  E,  D  (37) 

■*  °i 

where  the  colored  noise  is  determined  by: 


h  s  -  r-  5i  +  V 

Si  *1 


(38) 


The  t  represents  the  correlation  time  of  the  colored  noise 
®i 

and  the  v  the  strength  of  the  driving  white  noise  obtained 
^  2 

using  the  specified  variance  a..  of  the  colored  noise  and 

5  X 

the  formula 


*  *,1 'V 


(39) 


The  quantity  w  is  a  white  noise  of  specified  strength 


«i 
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2.  Accelerometer  Measurement  Uncertainties 


The  accelerometer  measurement  uncertainties  (u)f^ 
and  (u)fE  are  modeled  in  the  same  way  as  the  gyro  drift 
uncertainties,  as  colored  noise  plus  white  noise: 

(u)f.  -a.  +  w  ;  i  =  N,  E  (MO) 

1  i  ai 

where  again  w  is  the  wnite  noise  of  specified  strength 
ai 

and  the  colored  noise  is  given  by: 


*  -  TT  at  *  vat 


(Ml) 


where  t  is  the  correlation  time  of  the  colored  noise 

with  variance  j  ,  and  the  strength  of  the  driving  white 
ai 

noise  v_  is  given  by: 
ai 


\  ‘  2  *•!  /t«l 


(M2) 


3*  Computer  Simulation  Results 

Using  the  same  set  of  equations  (35-a)  through 
(35~g)  but  introducing  the  appropriate  state  augmentation 
in  order  to  incorporate  the  exponentially  correlated  noise 
for  the  gyro  drift  and  the  accelerometer  measurement,  we 
simulated  the  operation  of  the  combined  I.  N.  S./G. P. S.  system 
and  achieved  the  following  result. 


The  operation  of  the  system  proved  to  be  excellent 
for  all  the  used  correlation  times  from  60  seconds  up  to 
3600  seconds  (1  hour).  The  attitude  and  navigation  errors 
were  found  to  behave  in  the  same  way  being  minimized  after  a 
period  of  one  hour.  Furthermore,  the  variation  of  the 
attitude  and  navigation  errors  is  similar  with  the  case  of 
the  white  noise  driven  I. N. S./G. P. S.  combined  system  which 
again  is  similar,  if  not  exactly  the  same,  with  the  ideal 
I. N. S./G. P. S.  system. 

In  Figures  94  through  100  we  present  the  I.  N.  S./G. P. S. 
system  operation  for  an  exponentially  correlated  input  noise 
with  a  correlation  time  of  1  hour  (3600  sec).  We  can  easily 
see  in  these  figures  that  the  behavior  of  the  combined 
I. N.S./G. P.S.  system  is  the  same  witn  that  when  white  noise 
drives  the  input  except  for  a  very  small  and  negligible 
increase  of  the  attitude  and  navigation  errors  after  the 
first  hour  of  operation. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

From  knowledge  gained  throughout  this  work  and  based  on 
the  material  presented  in  our  study,  the  following 
conclusions  are  drawn: 

As  far  as  the  I.N.3.  errors  are  concerned  we  saw  that 

1.  The  effects  of  constant  gyro  drift  errors  for  the 
stationary  case  are  related  to  Foucault  modulation  which  has 
only  a  second-order  effect  on  the  longitude  and  latitude 
error  states  and  permit  us  to  neglect  it  in  cheap  systems 
designed  for  navigational  purposes. 

2.  For  the  case  of  easterly  flight,  latitude  errors 
were  reduced  by  a  factor  of  1.557  which  correponds  to  the 
ratio  x/u>ie  and  Foucault  modulation  period  reduced 
analogously  from  33*9  hours  to  21.8  hours. 

3.  For  the  westerly  flight,  case  longitude  and  latitude 
errors  grow  in  approximate  proportion  to  the  time-drift  rate 
prod uct . 

The  accelerometer  bias  errors  have  the  same  effects 
for  the  stationary  case  as  the  gyro  drift  errors. 

5.  The  Foucault  modulation  period  increased  by  the  same 
as  above  factor  of  1.557  for  the  easterly  flight  case,  while 
for  the  westerly  flight  decreased  by  a  factor  of  0.442 
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corresponding  again  to  the  ration  x/u>ie  for  the  specific 
case  . 

6.  The  error  sensitivities  remain  unchanged  for  the 
easterly  and  westerly  flight  and  again  we  may  neglect 
Foucault  modulation  as  producing  only  second-order  effects 
on  the  navigation  states. 

For  the  combined  I.  N.  S. /G . P. S.  system  the  results 
achieved  by  this  study  proved  that  the  errors  of  the 
system's  state  variables  are  damped  out  in  less  than  one 
hour,  denoting  effective  and  successful  operation  of  the 
G.P.S.  aiding  to  the  I.N.S.  Specifically: 

7.  Using  suboptimal  Kalman  filter  gains  for  one  hour 
process,  the  Y-position  error  reduced  from  its  initial  value 
by  a  factor  of  6  in  one  hour  and  by  a  factor  of  17  in  four 
hour  s  . 

8.  With  the  same  suboptimal  Kalman  filter  gains  of  one 
hour  process,  the  X-position  error  proved  the  system  more 
attractive  since  the  error  reduced  by  a  factor  of  8.5  after 
one  hour  and  by  a  factor  of  856  after  four  hours. 

9«  Both  the  X  and  Y-velocity  errors  damped  out  very 
quickly  so  that  after  one  hour  the  Y-velocity  error  reduced 
by  a  factor  of  312.5  and  the  X-velocity  error  by  a  factor  of 
117.6,  while  for  a  four-hour  process  both  velocity  errors 
reduced  by  a  factor  of  571  from  its  initial  value. 
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10.  The  consideration  of  long  term  filter's  operation 
proved  no  divergence  at  all  for  a  process  of  36  whole  hours. 
The  errors  remained  at  the  same  attractive  levels  as  for  the 
four-hour  process,  fact  which  enables  us  to  conclude  that 
the  combined  I. N.S./G. P. S.  system  works  with  excellent 
results  for  both  sho'-t  and  long  term  periods. 

11.  Finally  the  operation  of  the  combined  I . N. S . /G. ?. S . 
system  under  exponentially  correlated  input  noise  proved  to 
be  excellent  for  all  different  correlation  times  from  60  sec 
up  to  3600  sec,  with  a  negligible  increase  in  the  attitude 
and  navigation  error  magnitudes  after  the  first  hour  of 
operation  . 


3.  SECOMMENDATIOMS 

Continued  study  of  this  work  can  be  based  on  the 
following  recommendations: 

1.  A  Kalman  filter  design  study  where  the  primary 
emphasis  will  be  placed  upon  determination  of  the  "best" 
filter  state  variable  vector.  A  general  covariance  analysis 
program  for  the  analysis,  evaluation  and  design  of  Kalman 
filters,  which  will  help  this  study,  has  been  tape  recorded 
from  the  Wright  Patterson  Air  Force  Base,  Air  Force  Avionics 
Laboratory  and  modified  by  the  author  for  use  in  N.P.S. 
campus  computer  . 

2.  Investigation  of  various  measurement  rates  using  the 
external  range  measurements  from  a  set  of  satellites  in  view 
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among  the  18  of  the  G.P.S.  and  Kalman  filter's  performance 
for  these  rates. 

3.  Possible  use  of  a  flight  profile  generator  program, 
which  will  generate  simulated  flight  patterns  instead  of 
considering  specific  only  cases  for  stationary,  easterly, 
and  westerly  flights,  together  with  a  satellite  motion 
generator  required  to  provide  necessary  information  regard¬ 
ing  the  satellites'  orbital  elements.  This  recommendation 
applies  only  to  U.S.  citizens  since  such  programs  already 
exist  but  they  are  classified. 

4.  Investigation  and  results  evaluation  for  the  effects 
upon  filter  performance  when  range-rate  measurements  are 
available.  Then  a  comparison  with  the  usage  of  only  range 
measurements  could  be  extracted.  Another  aspect  for 
investigation  could  be  the  satellite  bearing  measurements  to 
declare  best  observable  satellites  and  to  provide  better 
accuracy. 

5.  Finally  a  comparison  of  sequential  versus  simul¬ 
taneous  measurement  would  be  another  area  of  interest.  The 
performance  of  a  filter  working  with  sequential  measurements 
is  of  interest  primarily,  because  of  the  increased  cost  of 
equipment  -  -nuired  to  perform  simultaneous  measurements  and 
computations  as  compared  to  the  sequential  ones. 
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X-SCRLE=  1  .  QOE  +  O  1  UNITS  INCH,  [hours] 

T-SCRLE  =  5. 00E-Q5  UNITS  INCH,  [rad/meru] 

KWSTRS 

RUN  1  E  'N  '  V5  TIME 

Figure  4.  Stationary  Case.  North  Level  Error  [rad/meru] 
for  Constant  North  Gyro  Drift  [1  meru] . 
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X-SCRLE=1  .00E+01  UNITS  INCH,  [hours] 

Y -SCALE  =  5 . 00E-05  UNITS  INCH.  [Rad/meru] 

KNSTflS 

RUN  1  E 'E 1  VS  TIME 


Figure  5.  Stationary  Case.  East  Level  Error  [Rad/meru]  for 
Constant  North  Gyro  Drift  [1  meru] . 
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X  —  SCRlE  —  1  .  OOE  +  O  1  UNITS  INCH,  [hours] 

T-SCflLE  =  S.  00E-U4  UNITS  INCH.  [ Rad/me ru] 

KWSTRS 

RUN  1  E 'D ’  V5  T IME 

Figure  6.  Stationary  Case.  Azimuth  Level  Error  [Rad/meru] 
for  Constant  North  Gyro  Drift  [1  meru] . 
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X  -SCALE=  1  .  OOE  +  0  1  UNITS  INCH,  [hours] 
T-SCflLE=5. 00E-04  UNITS  INCH.  [Rad/meru] 

KWSTflS 

RUN  1  DLh  VS  TINE 

Figure  7.  Stationary  Case.  Latitude  Error  [Rad/meru]  for 
Constant  North  Gyro  Drift  [1  meru] . 
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X-SCRLE= 1 . OOE  +  O  1  UNITS  INCH,  [hours] 
T-SCRLE=2. 00E-03  UNITS  INCH.  (Rad/meru] 

KWSTflS 

RUN  2  DLO  V 5  TIME 


Figure  8.  Stationary  Case.  Longitude  Error  [Rad/meru]  for 
Constant  North  Gyro  Drift  [1  meru] . 
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X-SCRLE=1.00E+01  UNITS  INCH,  [hours] 

Y -SCALE  =  5 . 0QE-Q4  UNITS  INCH.  [Rad/hour «meru] 

KWSTRS 

RUN  2  DLOD  VS  TIME 


X-SCfllE=  1  .  OOE  +  O  1  UNITS  INCH,  [hours] 

T  -SCALE  =  1 . 00E  -0  1  UNITS  INCH.  [Rad/meru] 

KWSTfiS 

RUN  1  E ' N ’  VS  TIME 

Figure  11.  Easterly  Flight  at  600  ft/sec. 

North  Level  Error  [Rad/meru]  for  Constant 
North  Gyro  Drift  [1  meru] . 
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X-SCRLE=  1  .  00E  +  0  1  UNITS  INCH,  [hours] 

Y  -SCRLE=  1  .  OOE-O  1  UNITS  INCH.  [Rad/meru] 

KWSTflS 

RUN  1  E  'E  '  VS  T 

Figure  12.  Easterly  Flight  at  600  ft/sec 

East  Level  Error  [Rad/meru]  for  Constant  North 
Gyro  Drift  [1  meru]. 
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X-3CRLE=1.00E+QI  UNITS  INCH,  [hours] 
Y-SCRLE=2.00E+Q0  UNITS  INCH.  [Rad/meru] 

KWSTflS 

RUN  1  E  ’  0  ’  VS  TIME 

Figure  13.  Easterly  Flight  at  600  ft/sec. 

Azimuth  Level  Error  [Rad/meru]  for  Constant 
North  Gyro  Drift  [1  meru] . 
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i 


X -SCPLE= 1 . OQE  +  O  1 
T-SCflLE=l . OOE+OO 


K  N  5  T  fl  S 
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Figure  14.  Easterly  Flight  at  600  ft/sec. 

Latitude  Error  [Rad/meru]  for  Constant  North 
Gyro  Drift  [1  meru]. 


ME 


94 


I 


X -SCALE= 1 . OOE  +  O  1 
T-SCALE=1.00E+01 


KWSTflS 
RUN  2 


UNITS 

UNITS 


INCH,  [hours  ] 

INCH.  [Rad/meru] 

D L 0  VS  TIME 


Figure  15.  Easterly  Flight  at  600  ft/sec. 

Longitude  Error  [Rad/meru]  for  Constant  North 
Gyro  Drift  [1  meru]. 
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Figure  16.  Easterly  Flight  at  600  £t/sec. 

Latitude  Rate  Error  [Rad/hour*meru]  for  Constant 
North  Gyro  Drift  [1  meru]. 
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Figure  17.  Easterly  Flight  at  600  ft/sec. 

Longitude  rate  error  [Rad/hour • meru]  for  Constant 
North  Gyro  Drift  {1  meru] . 
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Figure  18.  Westerly  Flight  at  600  ft/sec. 

North  Level  Error  [Rad/meru]  for  Constant 
North  Gyro  Drift  [1  meru] . 
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Figure  20.  Westerly  Flight  at  600  ft/sec. 

Azimuth  Level  Error  [Rad/meru]  for  Constant 
Notth  Gyro  Drift  [1  meru]. 
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Figure  21.  Westerly  Flight  at  600  ft/sec. 

Latitude  Error  [Rad/meru]  for  Constant  North 
Gyro  Drift  [1  meru] . 


101 


X -SCRLE=  i  .  OQE  +  Q  1  UNIT'S  INCH. 
T-SCflLE=5.00E+QG  UNITS  INCH. 

KWSTflS 

RIJN  2  DLL 


[hours] 

[Rad/meru] 


TIME 


Figure  22.  Westerly  Flight  at  600  £t/sec. 

Longitude  Error  [Rad/meru]  for  Constant  North 
Gyro  Drift  [1  mem]. 
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Figure  23.  Westerly  Fligltat  600  ft/sec. 

Latitude  Rate  Error  [Rad/hour • me ru]  for  Constant 
North  Gyro  Drift  [1  meru] . 
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Figure  24.  Westerly  Flight  at  600  ft/sec. 

Longitude  Rate  Error  [Rad/hour*meru]  for  Constant 
North  Gyro  Drift  [1  meru]. 
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Figure  25.  Stationary  Case.  North  Level  Error  [Rad/200Ug] 
for  Constant  North  Accelerometer  Bias  [200yg] . 
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Figure  26.  Stationary  Case.  East  Level  Error  [Rad/200pg] 
for  Constant  North  Accelerometer  Bias  [200yg]. 
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Figure  27.  Stationary  Case.  Azimuth  Level  Error 

[Rad/2Q0yg]  for  Constant  North  Accelerometer 
Bias  [200yg] . 


107 


1 


Ou)  001  002  003  0C4 
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RUN  1  DLfl  VS  TIMS 

Figure  28.  Stationary  Case.  Latitude  Error  [Rad/200yg]  for 
Constant  North  Accelerometer  Bias  [200vg]  . 
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X-SChLE=  1  .  QuE+O  1  UNITS  INCH,  [hours] 
T-SCRLE=  1  .  00E-06  UNITS  INCH.  [Rad/200yg] 

K  W  5  T  fl  S 

RUN  2  DLO  VS  TIME 

Figure  29.  Stationary  Case.  Longitude  Error  [Rad/200yg] 
for  Constant  North.  Accelerometer  Bias  [200ygJ  . 
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Figure  30.  Stationary  Case.  Latitude  Rate  Error 

[Rad/hour* 200yg]  for  Constant  North  Accelerometer 
Bias  [2 OOpg ] . 
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KWSTflS 

RUN  2  DLOD  VS  TIME 

Figure  31.  Stationary  Case.  Longitude  Rate  Error 

[Rad/hour* 200yg]  for  Constant  North  Accelerometer 
Bias  [200yg] . 
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Figure  35.  Easterly  Flight  at  600  ft/sec. 

Latitude  Error  [Rad/200yg]  for  Constant 
North  Accelerometer  Bias  [200yg] . 
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Y -SCALE  =  1 . 00E -OS  UNITS  INCH.  [Rad/200Pg] 

KNS  TflS 

RUN  2  DLO  VS  T 

Figure  36.  Easterly  Flight  at  600  ft/sec. 

Longitude  Error  [Rad/200wg]  for  Constant  North 
Accelerometer  Bias  [200yg]  . 
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Figure  37.  Easterly  Flight  at  600  ft/sec. 

Latitude  Rate  Error  [Rad/hour • 200yg]  for  Constant 
North  Accelerometer  Rias  f200pg]. 
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Figure  38.  Easterly  Flight  at  600  ft/sec. 

Longitude  Rate  Error  [Rad/hour*200ug]  for  Constant 
North  Accelerometer  Bias  f200vg]. 
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Figure  39.  Westerly  Flight  at  600  ft/sec. 

North  Level  Error  [Rad/200ug]  for  Constant 
North  Accelerometer  Bias  [200pg] . 
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Figure  40.  Westerly  Flight  at  600  ft/sec. 

East  Level  Error  [Rad/200pg]  for  Constant  North 
Accelerometer  Bias  [200yg]. 
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X -SC  ALt  = 1 . OQE  +  G 1  UNITS  INCH,  [hours] 
Y-SCALE=5. 00E-07  UNITS  INCH.  [Rad/200ng] 

KWSTflS 

RUN  1  E  ’  □  ’  V S 

Figure  41.  Westerly  Flight  at  600  ft/sec. 

Azimuth  Level  Error  [Rad/200yg]  for  Constant 
North  Accelerometer  Bias  [200pg] . 
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Figure  42. 


Westerly  Flight  at  600  ft/sec. 

Latitude  Error  [Rad/200yg]  for  Constant  North 
Accelerometer  Bias  [200yg] . 
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X-SCALE=  i  .  OQE  +  Q  1  UNITS  INCH,  [hours] 

Y -5CALE= 1 . CGE-QS  UNITS  INCH.  [Rad/200ug] 

KWSTRS 
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Figure  43.  Westerly  Flight  at  600  ft/sec. 

Longitude  Ettot  [Rad/200yg]  for  Constant  North 
Accelerometer  Bias  [200ng], 
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Y-3CRLE  =  2. 00E-C6  UNITS  INCH  .  [Rad/hour  •  200  m  g] 

KW5TR5 

RUN  2  DLfiD  V5  TIME 

Figure  44.  Westerly  Flight  at  600  ft/sec. 

Latitude  Rate  Error  [Rad/hour*200ug]  for  Constant 
North  Accelerometer  Bias  [200yg]. 
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Y-SCRLE  =  5  .  OQE-OS  UNITS  INCH.  [Rad/houx ,  200  ugj 
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RUM  2  DLOD  VS  TIM 


Figure  45.  Westerly  Flight  at  600  ft/sec. 

^nfitXide  ErTOr  tRad/h°ur.  200yg]  for  Constar 

North  Accelerometer  Bias  [200vg]. 
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Figure  46.  Stationary  Case.  North  Level  Error 

[Rad/140yradJ  for  Initial  North  Level  Error 
[14Gyrad] . 
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Figure  48.  Stationary  Case.  North  Level  Error  [Rad/140yrad] 
for  Initial  East  Level  Error  [140prad]. 
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Y -SCfiLE =5 . COE -02  UNITS  INCH.  [Rad/i40prad] 
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Figure  49. 


Stationary  Case.  East 
for  Initial  East  Level 


Level  Error  [Rad/140yrad] 
Error  [140yrad] . 
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KWSTflS 

RUN  5  E  '  N  '  VS  TIN 

Figure  50.  Stationary  Case.  North  Level  Error  [Rad/140yradJ 
for  Initial  Azimuth  Level  Error  [140yrad] . 
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Figure  51.  Stationary  Case.  East  Level  Error  [Rad/140yrad] 
for  Initial  Azimuth  Level  Error  [140yrad] . 
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KWSTflS 
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Figure  52.  Stationary  Case.  North.  Level  Error 

[Rad/ (2  ft/sec)]  for  Initial  Latitude  Rate  Error 
[0 . 345nrad/hour  =  7  ft/sec]. 
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Figure  53.  Stationary  Case.  Azimuth  Level  Error 

[Rad/ (2  ft/sec)]  for  Initial  Latitude  Rate  Error 
[0.345  mrad/hour  *  2  ft/sec]. 
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Figure  54.  Stationary  Case.  Latit  ide  Error 

[Rad/(2  ft/secl]  for  Initial  Latitude  Rate  Error 
[0 . 345mrad/hour  *  2  ft/sec]. 
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Figure  55.  Stationary  Case.  Longitude  Error 

[Rad/(2  ft/sec) ]  for  Initial  Latitude  Rate  Error 
[0 . 345rarad/hour  a  2  ft/sec]. 
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Figure  56.  Stationary  Case.  Latitude  Rate  Error 

[ (Rad/hour) / (2  ft/sec)]  for  Initial  Latitude 
Rate  Error  [0.345  mrad/hour  »  2  ft/sec]. 
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Figure  57. 


Stationary  Case.  Longitude  Rate  Error 
[(Rad/hour)/ (2  ft/sec)]  for  Initial  Latitude 
Rate  Error  [0 . 345mrad/hour  *  2  ft/sec]. 
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Figure  58.  Stationary  Case.  North  Level  Error 

[Rad/(2  ft/sec) ]  for  Initial  Longitude  Rate  Error 
[0. 345mrad/hour  3  2  ft/sec] . 
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Figure  59.  Stationary  Case.  East  Level  Error 

[Rad/(2  ft/sec)]  for  Initial  Longitude  Rate 
Error  [0 . 345mrad/hour  =  2  ft/sec]. 
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Figure  60.  Stationary  Case.  Latitude  Error  [Rad/ (2  ft/sec)] 
for  Initial  Longitude  Rate  Error 
[0 . 34 5m rad/hour  *  2  ft/sec]. 
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Figure  61. 


Stationary  Case.  Longitude  error  [Rad/[2  ft/sec'l 
for  Initial  Longitude  Rate  Error  7  CJ 
[0 . 34Smrad/hour  ~  2  ft /sec]. 
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Figure  62.  Stationary  Case.  Latitude  Rate  Error 

[ ( Rad/hour) / (2  ft/sec)]  for  Initial  Longitude 
Rate  Error  [0 . 345mrad/hour  *  2  ft/sec]. 
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Figure  63.  Stationary  Case.  Longitude  Rate  Error 

[(Rad/hour) / (2  ft/sec)]  for  Initial  Longitude 
Rate  Error  [0 . 345mrad/hour  a  2  ft/sec]. 
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Figure  65.  Theoretical  I .N. S . /G.P.S. ,  4-Hour  Process. 

Fast  ’tritude  Error  [Rad]. 
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Figure  66.  Theoretical  I .N.S./G.P.S. ,  4-Hour  Proces 
Azimuth  Attitude  Error  [Rad]. 
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Figure  67.  Theoretical  I.N.S./G.P.S. ,  4-Hour  Process. 

v  Position  Error  [Rad]. 
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Figure  68.  Theoretical  I. N. S./G. P. S. ,  4-Hour  Proce; 
X  Position  Error  [Rad]. 
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Figure  69.  Theoretical  I.N.S./G.P.S. ,  4-Hour  Process. 
Y  Velocity  Error  [Rad/hour] . 
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Figure  70.  Theoretical  I .N.S. /G.P. S. ,  4-Hour  Process. 

X  Velocity  Error  [Rad/hour  1 . 
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Figure  71.  Normalized  Input  Noise  Versus  Time. 
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Figure  72.  Normalized  Measurement  Noise  versus  Time, 
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Figure  73.  Theoretical 
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P.S.,  36-Hour  Process. 
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Figure  76.  Theoretical  I.N.S./G.P.S. ,  56-Hour  Process. 

Y  Position  Error  [Rad]. 
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Figure  77.  Theoretical  I.N.S./G.P.S. ,  36-Hour  Process. 
X  Position  Error  [Rad]. 
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Figure  73.  Theoretical  I.N.S./G.P.S. ,  36-Hour  Process. 
Y  Velocity  Error  [Rad/hour]. 
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Figure  79.  Theoretical  I.N.S./G.P.S. ,  36-Hour  Process. 

X  Velocity  Error  [Rad/hour] . 
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Figure  80.  Realistic  I.N.S./G.P. S. ,  4-Hour  Process. 
North  Attitude  Error  [Rad] . 
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Figure  81.  Realistic  I.N.S./G.P.S. ,  4-Hour  Proces 
East  Attitude  Error  [Rad] . 
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Figure  82.  Realistic  I .N.S ./G. P.S. ,  4-Hour  Process. 

Azimuth  Attitude  Error  [Rad] . 


162 


..£J0 


•  icmle  =  :  . ::e-co  jni’s  inch,  [hours] 

:  '  -:'_E  =  5  .  '  nr. -  C5  U N  I  T S  I\CU.  [Rad] 

•'  \  3  x  --  s 

’j\  :  j_-T  Vb  T 

Figure  83.  Realistic  I.N.S./G.P.S. ,  4-Hour  Process. 

Y  Position  Error  [Rad]. 
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Figure  84.  Realistic  I.N.S./G.P.S. ,  4-Hour  Process. 
X  Position  Error  [Rad], 
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Figure  85.  Realistic  I.N.S./G.P.S. ,  4-Hour  Process 
Y  Velocity  Error  [Rad/hour] . 
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Figure  86.  Realistic  I.N.S./G.P.S.  ,  4-Hour  Process. 

X  Velocity  Error  [Rad/hour] . 
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Figure  87.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process 
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Figure  88.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process. 
East  Attitude  Error  fRad]  . 
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Figure  89.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process. 
Azimuth  Attitude  Error  [Rad] . 
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Figure  90.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process. 
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Figure  91.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process 
X  Position  Error  [Rad]. 
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Figure  92.  Realistic  I.N.S./G. P.S. ,  36-Hour  Process. 

Y  Velocity  Error  [Rad/hour]. 
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Figure  93.  Realistic  I.N.S./G.P.S. ,  36-Hour  Process. 

X  Velocity  Error  [Rad/hour] . 
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Figure  94.  Realistic  I.N.S./G.P.S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  North  Attitude 
Error  [Rad] . 
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Figure  95.  Realistic  I.N.S./G.P.S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  East  Attitude 
Ertor  r  Rh  d 1 
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Figure  96.  Realistic  I.N.S./G.P.S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  Azimuth  Attitude 

Error  [Rad]. 
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Figure  97.  Realistic  I .N. S./G.P.S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  Y-Position  Error 
[Rad] . 


X-SCRLE=  !  .  OOE  +  OO  UNITS  INCH,  [hours] 

Y-SCALE=  1  .  00E-04  UNITS  INCH.  [Rad] 

KWSTflS 

RUN  2  DLQNG  VS  TIME 

Figure  98.  Realistic  I .N. S. /G. P. S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  X-Position  Error 
[Rad] . 
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Figure  99.  Realistic  I.N. S./G. P. S. ,  Exponentially  Correlated 
Input  Noise.  4 -Hour  Process.  Y-Velocity  Error 
[Rad/h.our] . 
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Figure  100.  Realistic  I .N. S. /G. P. S. ,  Exponentially  Correlated 
Input  Noise.  4-Hour  Process.  X-Velocity  Error 
[Rad/hour] . 
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’  V. 


APPENDIX  A 


A  SIMPLE  EXAMPLE: 

KALMAN  FILTER  APPLICATION  TO  A  RADAR  POSITION 
AIDED  INERTIAL  NAVIGATION  SYSTEM 

I.  INTRODUCTION 

The  application  of  3  Kalman  filter  to  3  simplified  radar 
position  aided  inertial  navigation  system  was  investigated 
as  a  first  step  of  our  study.  Since  the  case  appears  to  be 
easy  to  understand  difficult  concepts  and  the  results  prove 
the  design  expectations  we  present  hereafter  this  simple 
case  formulated  according  to  the  concepts  and  the  outline  of 
Ref.  2. 


The  I.N.3.  system  was  modeled  as  white  noise  driving  a 
0 

1/s  plant.  Radar  measurements  were  assumed  to  be  availaole 
and  were  similarly  corrupted  by  white  noise. 

The  differential  equations  describing  the  system  and  the 
Kalman  filter  were  numerically  integrated  to  yield  the 
response  for  a  wide  range  of  input  conditions  and  system 
noi3e  statistics.  Particular  attention  was  paid  to 
conditions  in  which  the  noise  statistics  employed  in  the 
filter  were  different  from  the  statistics  of  the  noise 
actually  driving  the  system  dynamics  and  measurements. 

For  all  cases  considered,  including  those  for  which 
intentional  mismatches  in  the  noise  statistics  were 
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introduced,  the  filter  was  found  to  perform  in  an  entirely 
satisfactory  manner.  This  is  the  filter  reliably  and  quite 
accurately  tracked  the  system's  dynamics  even  at  the 
presence  of  at  times  rather  severe  levels  of  noise. 

In  the  section  to  follow,  the  theoretical  development 
of  the  Kalman  filter  equations  will  be  presented.  This 
development  is  based  on  that  given  in  Chapter  6  of  Maybeck 
[Ref. 2]  and  according  to  explanations  given  in  class  notes 
from  Prof.  Collins  [Ref.  12]. 

Next,  a  discussion  of  simulation  results  will  be  given, 
in  which  the  various  cases  considered  are  outlined,  and  the 
performance  of  tne  filter  in  each  case  is  described. 

Finally  an  overall  summary  and  conclusions  regarding  the 
oosarved  performance  of  the  filter  over  a  wide  range  of  test 
conditions,  is  presented. 

II.  I.N.S.  AIDED  BY  POSITION  DATA 

For  this  problem,  the  model  of  the  I.N.S.  is  simply  a 
double  integration  of  noise-corrupted  acceleration  infor¬ 
mation,  as  depicted  in  Fig.  101.  The  noise  w  is  a  white 
Gaussian  noise  of  zero  mean  and  variance  Kernel 

£[w( t)w( t+T>]  =  Q«(t) 

entering  at  the  acceleration  level,  and  it  is  meant  to  model 
the  errors  corrupting  the  I.N.S.  accelerometer  outputs 
(accelerometer  biases  and  noise,  platform  misalignment, 
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etc.)*  The  noise-corrupted  acceleration  is  integrated  once 
to  yield  I. N. S. -indicated  velocity  (v^),  and  a  second  time 
to  obtain  inertially-indlcated  position  ( r i ) . 

Similarly  a  simple  model  for  the  radar  or  radio 
navigation  aid  is  the  true  position  (r^)  corrupted  by  noise 
u,  which  is  again  white  Gaussian  with  zero  mean. 

The  two  error  state  variables  for  this  case  are: 

sr  ( t)  =  r .  ( t)  -  rt  ( t) 

1  L  (A-1) 

Sv(t)  =  V  i  ( t )  -  Vt(t) 

The  measurement  to  be  presented  to  the  filter  is  the 
difference  between  the  inertially  indicated  position  and 
that  measured  by  the  radar  or  radio  navigation  aid. 

From  Figure  3b  we  have: 

z(t)  =  r^t)  -  rr(t) 

=  l>t(t)  +  Sr  ( t)  ]  -  [rt(t)  -  u(t)]  =  (A-2) 

=  Sr ( t)  +  u(  t) 

Thi3  is  a  measurement  of  the  error  Sr(t)  corrupted  by  noise 
u(  t)  . 

To  establish  the  state  dynamics  model  for  the  error 
states,  first  let  us  consider  the  total  states  r^  and  v^. 
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(6-3) 


The  true  position  velocity  and  acceleration  are  related 


by  : 


r  •  i 

r 

r  -1 

0  1 

rt 

0 

• 

V. 

0  0 

1 

U 

,  j 

L  -J 

-  - 

( A-4) 


Subtracting  (A-4)  from  ( A — 3 )  and  using  the  error  state 
definitions  of  (A-1)  we  find  the  desired  relation  as: 


• 

- 

*1 

r 

Sr 

0 

1 

Sr 

0 

5v 

0 

0 

57 

1 

.  « 

- 

■  J 

L  -1 

-  - 

(A-5) 


The  measurement  model  z  can  be  expressed  in  terms  of 
errors  as  : 


z(t)  =  [1  0] 


Sr 

5v 


+  u(t)  =  H  x  +  u  ( t ) 


Since  we  would  like  to  design  a  Kalman  filter  for  this 
situation  we  need  to  solve  the  RICATI  equation  as  below: 


P  =  FP  +  PFT  +  GQGT  -  PHTR“1HP 


( A-6 ) 


where 


P  = 


t 

1 1 

P  12 

21 

P22 

and 


P 12  =  P21 
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Kalman  Filter  Block  Diagram 


E[ w( t) W( t+T) ]  s  Q 5 ( T ) 
E[u(t)u(t+T)3  =  R«(t) 


Solving  for  RICATI  equation  for  our  case  we  get: 


P  P 

r11  12 


P  P 

21  22 


2P12  -  (Pn)VR 


P22  -  P11P12/R 


P22  '  P11P12/H 


Q  -  (P12)£/R 


(A-6a) 


For  the  steady  state  case  where  P  =  0  we  get  tne 

elements  of  covariance  matrix  in  terms  of  Q,  R  and  the  ratio 
1/4 


( Q/ R )  which  represents  the  natural  frequency  of  the 
system : 


Pn  -  2  Q 1  /4  R3/4 


12 


q1/2  ,1/2  s  p 


21 


(A-7) 


P22  =  2  Q3/4  R1/4 


The  Kalman  filter  equation  is 
x  :  Fx  +  Gu  +  K(t)  [z  -  H_x  ] 
which  in  terms  of  error  quantities  becomes: 


(A-8) 
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sr 

0  1 

sr 
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6  V 
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5  v 

+  K  [z  -  sr] 


( A-9 ) 
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where  K  =  PHTR’1 


‘p11 

P12 

1 

R21 

M. 

R22 

0 

m  - 

r-„/H 
p  1 2/  R 


1 

If 


(A-10) 


and  using  the  results  of  the  RICATI  equation  solution  we  can 
write 


1 

S2 

K2 

_  J 

(-n>2  _ 

where  [r*ad/sec]  ( A— 1 1 ) 


From  the  above  information  we  oan  draw  the  block  diagram  of 
the  Kalman  filter  for  this  system  as  shown  in  Figure  102. 

The  initial  transient  behavior  of  tne  filter  gains  < . 
and  K^  depends  on  PQ ,  but  they  are  within  a  few  percent  ot 
their  steady  state  values  (independent  of  PQ)  after  unt  =  2, 
so  a  prediction  of  time  to  reach  steady  state  would  be 
approximately  2/&>n  seconds  [Ref.  2]. 

The  filter  can  be  put  into  either  feedforward  configura¬ 
tion  or  feedback  configuration .  Since  for  our  study  we  use 
the  feedback  configuration  we  present  here  the  outline  and 
the  results  for  this  configuration .  A  block  diagram  of  the 
system  is  presented  in  Figure  103  as  depicted  in  Maybeck'3 
work  [Ref.  2].  This  block  diagram  allows  us  to  write  the 
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system’s  equations  which  we  will  simulate  numerically  to 
achieve  the  system's  performance. 

We  define  the  outputs  of  the  I.N.S.  system  corrected  by 
feedback  from  the  filter  as  follows: 


r(t)  =  ri(t)  -  «r(t) 
v ( t )  =  v.(t)  -  5  v  (  t ) 


( A-12) 


which  will  be  a  very  helpful  mathematical  tool  since  the 
most  straightfor ward  means  of  generating  feedback 
implementations  is  to  write  the  system  and  filter  equations 
in  terms  of  corrected  system  states.  (For  our  case 
corrected  I.N.S.  states.) 

Then  using  the  equations  ( A — 3 )  and  ( A-9 )  together  wicn 
equation  (A-8)  we  can  write  the  matrix  form  of  the  system 
equations  as: 


A 

rfc( t)-«r(t) 

■°  r 

r  j^CtJ-srlt) 

V 

[at ( t)+w(  tj 

-Kl(t)' 

A 

= 

- 

^(tl-Svft) 

0  0 

v^( t)-Sv( t) 

_1  _ 

K2(t) 

and  finally  we  get  the  simple  form: 


[  z(  t ) -Sr ( t )  ] 

CA-13) 


A 

r  ( t) 

- 

0  1 

r(t) 

0 

|at(  t)+w(  t)] 

K 1  ( t) 

[rp( t) -r ( t)  ] 

A 

s 

v(t) 

0  0 

_v(t)_ 

1 

K2(t) 

C  A— 14) 

In  the  next  section  the  programming  and  simulation 


results  of  the  system  i3  presented. 


gure  103.  Feedback  Kalman  Filter  Configuration 


III.  COMPUTER  SIMULATION  AND  RESULTS 

We  simulated  the  system  for  different  input  signals, 
different  noise  characteristics  (zero  mean  Gaussian  noise 
with  different  variances)  to  see  the  effect  of  the  filter 
for  error  estimation. 

From  the  block  diagram  of  the  system  in  Figure  103  we 
can  write  the  following  set  of  equations  which  we  will  use 
to  get  results  by  numerical  computer  simulation. 

x( 1 )  =  rfc  =  vt 

k{2)  =  vt  =  afc 

*<(3)  =  r.  =  v. 

X  (  4  )  =  V  .  =  3„  +  w 

x(5)  =  Pn  =  2P12  -  (Pn)2/R 

x(6)  =  P12  =  P22  -  Pnp22/R  ( A-1 5 ) 

x(7)  =  P22  =  Q  -  (P12)2/R 

A  A 

x(8)  =  Sr  s  Sv  +  K^(z  -  sr) 

x(9)  =  4v  =  K2(z  -  Sr) 

x( 15)  =  r p  =  rfc  -  u 

x  ( 16 )  =  v  =  v  ^  -  4  v 
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The  above  set  of  differential  equations  of  the  I.N.S. 
system  and  feedback  Kalman  filter  were  numerically 
integrated  using  INTZG2  computer  routine  in  conjunction  witn 
a  routine  (LNORM)  for  generating  Gaussian  distributed  random 
numbers  to  represent  the  noise  into  the  system.  Typical 
simulation  runs  used  an  integration  step  size  of  0.01 
seconds  and  a  total  run  time  of  36  seconds  by  which  point 
the  filter  had  easily  achieved  steady-state  operation  in 
most  cases. 

Shown  on  the  next  page  is  a  run  summary  representing  the 
various  conditions  that  were  tested.  For  each  of  six  cases, 
the  covariances  of  the  measurement  noise  (R)  and  process 
noise  (Q)  are  indicated.  Note  that  in  a  number  of  cases, 
the  noise  statistics  used  in  the  filter  were  chosen  to  be 
different  from  those  characterizing  the  input  noise  entered 
into  the  system.  This  intentional  "mismatch"  was  done  to 


investigate  filter  performance  under  conditions  where  the 
true  "real  world'*  noise  statistics  are  inadequately  or 
poorly  known . 

In  particular  it  was  desired  to  determine  whether  any 
instances  of  filter  "divergence"  could  be  observed  as  a 
result  of  the  mismatch  in  system  noise  statistics.  It  is 
noteworthy  that  for  all  conditions  tested,  the  filter 
performed  in  an  entirely  satisfactory  manner  with  no 
evidence  of  divergence. 

It  should  be  noted  here  that  in  Table  XI  the  noise 
covariances  Q  and  R  actually  represent  the  statistics  of 
Discrete  Noise  entering  the  system  at  the  integration 
interval  it  =  0.01  seconds.  That  is: 

Q  =  Etw^w.J] 

where  tk  =  kit 

R  =  ECu.uJ] 

As  it  is  pointed  out  by  Bryson  and  Ho  [Ref.  13]  the 
numerical  integration  is  a  discrete  approximation  to  a 
continuous  system  whose  noise  processes  have  spectral 
densities  given  by 

E[w(t)w(r+t) ]  =  Q ' 5 ( t ) 

E[ u( t) u( t+t) ]  =  R'«(t) 
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The  relation  between  Q  and  Q'  and  R  and  R'  is  according 
to  Bryson: 

Q*  =  Q*  At 

R'  =  R»  At 

Also  shown  in  Table  XI  for  each  case,  are  the  filter 
natural  frequency  and  the  steady-state  values  of  tne  Kalman 
gains . 

A  brief  discussion  will  now  be  given  of  the  results  for 
each  of  the  six  cases.  Detailed  plots  of  the  variables  of 
interest  for  each  case  are  attached  and  will  be  referred  to 
in  the  subsequent  discussion. 

1 .  Case  I 

for  this  case  we  used  R  =  10,000  and  Q  =  1  for 
the  filter.  The  actual  noise  however  is  mismatched  with 
=  and  Qt  3  ^  an<^  thus  the  filter  assumes  the 

measurement  noi3e  of  the  radar  position  data  considerably 
higher  than  the  case  is  actually.  Shown  in  the  attached 
plots  on  Figures  104  and  105  is  the  type  of  noise  actually 
entered  into  the  system  using  a  Gaussian  random  number 
generator.  Also  shown  are  the  histories  of  the  Kalman 
filter  gains  K-j  and  K2  versus  time.  The  performance  of  the 
filter  for  this  case  is  outstanding  as  evidenced  by  the  two 
plots  for  Case  I  in  Figures  108  and  109.  Here  the  estimated 
position  out  of  the  filter  coinsides  with  the  true  position 
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denoting  that  the  filter  tracks  the  system  extremely  well. 
Among  the  other  attached  plots,  Figure  112  presents  the 
noise  corrupted  radar  position  in  a  very  imposing  way. 

2.  Case  II 

This  case  represents  one  in  which  the  filter  and 
external  noise  are  "tuned”  so  that  the  same  noise  statistics 
are  employed  with  R  =  100  and  Q  =  25.  Again  the  Kalman 
gains  are  plotted  indicating  the  time  of  steady-state 
condition  in  Figure  114  and  115.  As  it  is  depicted  from 
Figures  116,  117  and  118  the  Kalman  filter  rapidly 
"locks-on"  to  the  true  position  and  velocity  and  accurately 
tracks  the  system  thereafter. 

3.  Case  III 

In  this  case  the  filter  and  external  noise  are 
"tuned"  with  R  =  100  and  Q  =  1.  The  system’s  initial 

p 

conditions  now  include  a  10  ft/sec  constant  acceleration 
and  it  was  desired  to  see  how  well  the  Kalman  filter  kept  up 
with  the  changing  input.  Once  again  the  performance  of  the 
filter  in  accurately  tracking  the  system  is  excellent.  This 
can  be  verified  looking  at  Figures  124  and  125  and  126  and 
127  respectively  where  we  can  see  the  coinsidence  of  the 
true  and  estimated  position  and  velocity. 

4.  Case  IV 

For  thi3  run  the  filter  and  external  noise  are  again 
"tuned"  but  with  increased  statistics  of  R  =  400  and  Q  =  50. 
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The  attached  set  of  plots  in  figures  1 3 1  through  138  present 
the  system  and  the  filter  operation  proving  the  accurate  and 
satisfactory  tracking  of  the  system 

5.  Case  V 

Now  the  filter  perceives  the  radar  measurement  noise 
to  be  higher  than  it  actually  is.  The  statistics  used  for 
this  case  were  R  =  400  and  Q  =  50  for  the  filter  while  for 
the  external  noise  we  used  R^  =  50  and  Qt  =  50.  The 
reliability  and  the  accuracy  of  the  system  is  again  depicted 
in  the  attached  plots  for  Case  V  in  Figures  139  through  146. 

6.  Case  VI 

In  the  last  case  considered  the  statistics  of  the 
random  process  noise  exciting  the  I.N.S.  accelerometers  was 
mismatched  witn  that  assumed  in  the  filter.  Here  we  used 
Qt  s  50  and  Q  =  10.  The  radar  measurement  noise  Rt  =  R  = 

400  was  assumed  the  same.  The  set  of  plots  in  Figures  147 
througn  154  indicate  very  good  performance  of  the  filter 
despite  the  intentional  mismatch  introduced  for  the  system 
driving  noise. 

IV.  CONCLUSIONS 

Following  the  development  of  Reference  2,  simplified 
equations  characterizing  the  Kalman  filter  were  derived  and 
numerically  integrated  to  yield  the  filter  response  to  a 
wide  range  of  input  conditions  and  system  noise  statistics. 
Particular  attention  was  paid  to  conditions  in  which  the 


noise  statistics  employed  in  the  filter  differed  from  the 
statistics  of  the  noise  records  actually  driving  the  system 
dynamics  and  measurements. 

For  all  cases  considered  including  those  for  which 
intentional  mismatches  in  the  noise  statistics  were 
introduced,  the  filter  was  found  to  perform  in  an  entirely 
satisfactory  and  reliable  manner.  By  that  is  meant  that  the 
filter  very  accurately  tracked  she  system  dynamics  even  in 
the  presence  of  at  times  rather  sever  levels  of  noise. 

Examination  of  typical  time  histories  for  the  variables 
of  interest,  showed  that  the  filter  Kalman  gains  and  K~ 
rapidly  settled  to  their  theoretical  steady  state  values 
within  a  time  short  compared  to  the  average  run  time.  This 
was  accompanied  oy  the  filter  range  and  velocity  estimates 
rapidly  locking  on  to  the  true  system  position  and  velocity 
and  accurately  tracking  it  thereafter. 

It  is  concluded  then  that  the  Kalman  filter  configura¬ 
tion  discussed  here  above  performed  extremely  well  over  the 
range  of  conditions  considered. 
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TABLE  VIII 


t 

'  COMPUTER  RUNS  SUMMARY 


Filter  True 


CASE 

R 

Cft]2 

Q 

sec 

Rt 

Cft]2 

0* 

sec 

u» 

,2  rad 
rj  T51T 

^SS 

(K2>ss 

I 

10, 000 

1 

100 

1 

0.  100 

0.  141 

0.010 

II 

100 

25 

100 

25 

0. 7  07 

1 

0.500 

III1 

100 

1 

100 

1 

0.  316 

0.447 

0.100 

IV 

400 

50 

400 

50 

0.595 

0.841 

0.354 

V 

400 

50 

50 

50 

0.595 

0.841 

0-354 

VI 

400 

10 

400 

50 

0.398 

0.562 

0.158 

Initi 

al  Cond 

itions  : 

Post  t ion 
Velocity 


S  200  ft 
=  50  ft/sec 


Acceleration  s  0 


1 1 


5  P 


22 


1A4 

=  10  , 


12 


21 


0 


R 

Q 

Rt 

Qt 

(K1>SS 

(K2)S3 


Natural  frequency  =  [Q/R]1^4 

Measurement  noise  covariance  used  in  filter 

Process  noise  covariance  used  in  filter 

True  measurement  noise  covariance 

True  process  noise  covariance 

Steady  3tate  gain  K1  =  2  u 

Steady  state  gain  K2  =  [«n32 


Vor  this  case  system  assumed 
acceleration  10  ft/sec  . 


to  have  constant 
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X-SCALE=  1 .  OOE  +  0 1  UNITS  INCH,  [hours] 

T-SCALE  =  2 . 00E  +  0  1  UNITS  INCH.  [ft/(sec)2] 

KWSTFIS 

RUN  1  N0ISE-W  VS  TIME 


Figure  104.  Normalized  Input  Noise  Versus  Time, 
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X-SCRLE=  1 .  OOE  +  0 1  UNITS  INCH,  [hours] 

T-SCRLE=2 . 00E  +  0 1  UNITS  INCH,  [(ft)2] 

KWSTflS 

RUN  1  NQISE-V  VS  TIME 

Figure  105.  Normalized  Measurement  Noise  Versus  Time. 
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X-SCRLE  =  1 . 00E+01  UNITS  INCH,  [hours] 
Y-SCPLE=5. 00E-01  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  1  K 1  VS 


Figure  106.  Case  I.  Kalman  Filter  Gain  to  Velocity. 


X-SCBLE=1.00E  +  01  UNITS  INCH,  [hours] 
Y-SCALE=2.00E-01  UNITS  INCH.  [ft/(sec)2] 

KWSTflS 

RUN  1  K2  VS  T I 


Figure  107.  Case  I.  Kalman  filter  gain  to  Acceleration 


X-SCRLE= 1 . OOE  +  O  1  UNITS 
Y-SCRLE=5. 00E+02  UNITS 


KWSTRS 
RUN  2 


INCH.  [hours] 

INCH,  [ft] 

R-TRUE  VS  TIME 


Figure  108.  Case  I.  True  Position  Versus  Time. 
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X-SCFtLE  =  1 . 00E  +  01  UNITS  INCH,  [hours] 
Y-SCFILE  =  5.0QE  +  02  UNITS  INCH,  [ft] 


KWSTflS 

RUN  2  R-HRT  VS  TIME 


Figure  109.  Case  I.  Predicted  Position  Versus  Time. 
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T 


X-SCRLE=1 . 00E  +  01  UNITS  INCH,  [hours] 

Y-SCRLE=2 . 00E  +  0 1  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  2  V-HRT  VS  TIME 


Figure  110.  Case  I.  Predicted  Velocity  Versus  Time, 
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X-SCRIE  =  1 .  OOE  +  O  1  UNITS  INCH,  [hours] 
Y-SCfllE=  t .  OOt  +  O  1  UNITS  INCH,  [ft] 


KWSTRS 

RUN  3  R-INS  VS  TIME 


Figure  111.  Case  I,  I.N.S.  Indicated  Position  Versus 


X-SCRLE=  1  .  OOE  +  O  1  UNITS  INCH,  [hours] 
Y-SCRLE=5.0QE+02  UNITS  INCH.  -  t] 

KWSTRS 

RUN  3  R-RRDRR  VS  TIME 

Figure  112.  Case  I,  Radar  Indicated  Position  Versus  Time 


X-SCRLE=1 . 00E+01  UNITS 
Y-$CRLE=2.00E+00  UNITS 


KWSTRS 
RUN  3 


INCH,  [hours] 

INCH,  [ft/sec] 

V-INS  VS  TIME 


Figure  113.  Case  I.  I.N.S.  Indicated  Velocity  Versus  Time 
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X-SCRLE=  1 . 00E  +  01  UNITS  INCH,  [hours] 
T-SCr)LE  =  2. 00E+0 1  UNITS  INCH.  [ft/Sec] 

KWSTflS 

RUN  1  K 1  VS 


Figure  114.  Case  II.  Kalman  Filter  Gain  to  Velocity. 


3  GOO  001  002  003  004 


X-SCRLE  =  1  .  OOE  +  Q  1  UNITS  INCH,  [hours] 

Y-SCALE  =  5.00E  +  00  UNITS  INCH.  [ft/(sec)2] 

KWSTflS 

RUN  1  K2  VS  TIME 


Figure  US.  Case  II.  Kalman  Filter  Gain  to  Acceleration. 


211 


Figure  116.  Case  II.  True  Position  Versus  Time. 


1 _ _ _ _ _ -■  -  : _  : 

000  OCi  002  003  00«4 


X-SCflLE=l . 00E+01  UNITS  INCH,  [hours] 
T-SCRLE=5.00E  +  0a  UNITS  INCH,  [ft] 


KWSTflS 

RUN  2  R-HflT  VS  TIME 


Figure  117.  Case  II.  Predicted  Position  Versus  Time. 
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X-SCRLE=  1  .  OOE  +  O  1  UNITS  INCH. [hours] 

T-SCRLE=2. Q0E+01  UNITS  I NCH .  [ft/sec] 

KWSTflS 

RUN  2  V-HRT  VS  TIME 

Figure  118.  Case  II.  Predicted  Velocity  Versus  Time. 


214 


I 


X-SCRLE*  1 . 00E  +  0 1  UNITS  INCH,  [hours! 

T-SCRLE  =  1.00E  +  01  UNITS  INCH,  [ft] 

KWSTflS 

RUN  3  R-INS  VS  TIME 


Figure  119.  Case  II.  I.N.S.  Indicated  Position  Versus  Time 


-SCfiLE=1.00E+01 

UNITS 

INCH. 

[hours] 

-SCRLE=5 . QQE+Q2 

UNITS 

INCH. 

[ft] 

KWSTRS 

RUN  3  R-RRDRR  VS  TIME 

Figure  120.  Case  II.  Radar  Indicated  Position  Versus  Time. 
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X-SCRLE=  1 .  OOE+0 1  UNITS  INCH,  [hours] 
Y-SCRLE=2.  OOE+OO  UNITS  INCH,  [ft/sec] 

KWSTRS 

RUN  3  V-INS  VS  TIME 

igure  121.  Case  II.  I.N.S.  Indicated  Velocity  Versus 


X-SCRLE= 1 . 00E+01  UNITS  INCH. 
Y-SCALE=2 .  00E+0 1  UNITS  INCH. 

KWSTRS 

RUN  1  K1 


[hours  ] 
[ft/sec] 

VS 


Figure  122.  Case  III.  Kalman  Filter  Gain  to  Velocity 


X-SCRLE  =  1 . 00E+01  UNITS  INCH,  [hours] 
T-SCRLE^S .  OOE+OO  UNITS  INCH.  [ft/(sec)z] 

KWSTRS 

RUN  1  K2  VS  T 


Figure  123.  Case  III.  Kalman  Filter  Gain  to  Acceleration 


X-SCRLE=1.Q0E  +  01  UNITS  INCH,  [hours] 

Y-SCRLE  =  5. 0QE  +  03  UNITS  INCH.  [£t] 

KWSTflS 

RUN  2  R-TRUE  VS  TIME 


Figure  124.  Case  III.  True  Position  Versus  Time. 
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000  005  0)0  015  020 


X-SCRLE=  1 .  OOE  +  O  1  UNITS  INCH,  [hoursj 
T-SCRLE=5.00E+03  UNITS  INCH,  [ft] 


KWSTRS 

RUN  2  R-HflT  VS  TIME 


Figure  125.  Case  III.  Predicted  Position  Versus  Time 


Figure  126.  Case  III.  True  Velocity  Versus  Time. 


X  —  SCFILE=1 .  OOE  +  Ol  UNITS  INCH,  [hours] 

T-SCRLE=1 . 00E+02  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  2  V-HRT  VS  TIME 


l 


Figure  127.  Case  III.  Predicted  Velocity  Versus  Time. 


000  002  004  006  008 


X-SCRLE  =  1 . 00E  +  01  UNITS  INCH,  [hours] 

Y-SCALE  =  2.00E  +  03  UNITS  INCH,  [ft] 

KWSTRS 

RUN  3  R-INS  VS  TIME 


Figure  128.  Case  III.  I.N.S.  Indicated  Position  Versus 
Time . 
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X-SCALE=  1 .  OOE  +  0 1  UNITS  INCH,  [hours] 
Y-SCfllE=5.00E+03  UNITS  INCH,  [ft] 

KWSTflS 

RUN  3  R-RRDRR  VS  TIME 


Figure  129.  Case  III.  Radar  Indicated  Position  Versus 
Time . 
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X-SCALE=1 . 00E+01  UNITS  INCH,  [hoursi 
T-SCALE  =  1 .  OOE+02  UNITS  INCH,  [ft/secj 

KWSTflS 

RUN  3  V-INS  VS  TIME 


Figure  130.  Case  III,  I.N.S.  Indicated  Velocity  Versus 

Time.  j 

I 

i 

i 
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X-SCRLE=  1 . 00E  +  01  UNITS  INCH,  [hours] 
T-SCHLE=5. 00E  +  00  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  1  K 1  VS 


Figure  151.  Case  IV.  Kalman  Filter  Gain  to  Velocity 


X-SCRLE=  1 .  OOE+O 1  UNITS  INCH.  [hours] 
T-SCRLE=2. 00E+00  UNITS  INCH.  [ft/(sec)2] 

KWSTflS 

RUN  1  K2  VS  T! 


Figure  132.  Case  IV.  Kalman  Filter  Gain  to  Acceleration 


X-SCflLE= 1 . QOE  +  O 1 
T-SCRLE=5.00E+02 


KWSTHS 
RUN  2 


UNITS  INCH.  [hours] 

UNITS  INCH,  [ft] 

R-TRUE  VS  TIME 


Figure  133.  Case  IV.  True  Position  Versus  Time. 
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X-SCfiLE  =  l .  00E  +  01  UNITS  INCH,  [hours] 

T-SCALE  =  5.  OOE  +  Oa  UNITS  INCH.  [ft) 

KWSTflS 

RUN  2  R-HRT  VS  TIME 


Figure  134.  Case  IV.  Predicted  Position  Versus  Time. 
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X-SCflLE=1.00E  +  01  UNITS  INCH.  [hours] 
T-SCfllE=2. 00E+01  UNITS  INCH.  [ft/sec] 

KWSTflS 

RUN  2  V-HRT  VS  TIME 


Figure  135.  Case  IV.  Predicted  Velocity  Versus  Time. 
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X-SCFILE  =  1 . 00E  +  01  UNITS  INCH,  [hours] 

Y -SCALE=2 . 00E  +  0 1  UNITS  INCH,  [ft] 

KWSTflS 

RUN  3  R-INS  VS  TIME 


Figure  136.  Case  IV.  I.N.S.  Indicated  Position  Versus 
Time . 


X-SChLE=  1 .  GOE  +  O  1  UNITS  INCH.  [hours] 
Y-SCfllE=5.00E+02  UNITS  INCH.  [ft] 


KWSTRS 

RUN  3  R-RAQRR  VS  TIME 


Figure  137.  Case  IV,  Radar  Indicated  Position  Versus 


X-SCflLE=1.00E  +  01  UNITS  INCH,  [hours] 

Y-SCflLE=5 . 00E+00  UNITS  INCH,  [ft/sec] 

KWSTRS 

RUN  3  V-INS  VS  TIME 

Figure  13.8.  Case  IV.  I.N.S.  Indicated  Velocity  Versus  Time. 
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X-SCHLE-1.00E+01  UNITS  INCH,  [hours] 
T-SCFILE  =  5.0QE  +  00  UNITS  INCH,  [ft/sec] 


KWSTHS 

RUN  1  K 1  VS 


Figure  139.  Case  V.  Kalman  Filter  Gain  to  Velocity. 


000  002  004  006  008 


X-SCFILE=1.00E  +  01  UNITS  INCH,  [hours] 
Y-SCRLE=2.00E+00  UNITS  INCH.  [ft/(sec)2] 

KWSTRS 

RUN  1  K2  VS  T 


Figure  140.  Case  V.  Kalman  Filter  Gain  to  Acceleration 


000  005  010  015  020 


I _ ; _ _ _ I _ I _ i _ 

000  001  002  003  004 


X-SCRLE  =  1 . 00E  +  01  UNITS  INCH,  [hours] 

Y-5CRLE=5.  QOE  +  02  UNITS  INCH,  [ft] 

KWSTflS 

RUN  2  R-TRUE  VS  TIME 


Figure  141.  Case  V.  True  Position  Versus  Time 


Figure  142.  Case  V.  Predicted  Position  Versus  Time. 


X-SCRLE=1 . 00E  +  01  UNITS  INCH,  [hours] 

Y-SCRLE=2. 00E+0 1  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  2  V-HflT  VS  TIME 


Figure  143.  Case  V.  Predicted  Velocity  Versus  Time. 


X-SCRLE= 1 . 00E+01 
Y-SCRLE=5.QQE+00 


KNSTRS 
RUN  3 


UNITS  INCH,  [hours] 

UNITS  INCH,  [ft] 

R-INS  VS  TIME 


Figure  144.  Case  V.  I.N.S.  Indicated  Position  Versus  Time. 
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X-SCflLE=l .  00E  +  01  UNITS  INCH,  [hours] 

T-SCflLE=5. 00E  +  02  UNITS  INCH,  [ft] 

KWSTflS 

RUN  3  R-RRDRR  VS  TIME 


Figure  145.  Case  V.  Radar  Indicated  Position  Versus  Time. 
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X-SCALE= 1  .  OOE+O 1  UNITS  INCH.Ihoursj 
T-SCALE=1  .  OOE+OO  UNITS  INCH,  [ft/sec] 

KWSTflS 

RUN  3  V-INS  VS  TIME 


Figure  146.  Case  V.  I.N.S.  Indicated  Velocity  Versus  Time. 


X-SCRLE=  1 .  OOE+O 1  UNITS  INCH,  [hours] 
T-SCALE  =  5.00E  +  00  UNITS  INCH,  [ft/sec] 

KWSTRS 

RUN  1  K 1  VS 


Figure  147.  Case  VI.  Kalman  Filter  Gain  to  Velocity 


X-SCRLE=1 . 00E+01  UNITS  INCH,  [hours] 
Y-SCALE=2.00E+00  UNITS  INCH.  (ft/(sec)2] 

KWSTRS 

RUN  1  K2  VS  T 


Figure  148.  Case  VI.  Kalman  Filter  Gain  to  Acceleration 


i 


X-SCALE=  1 . 00E  +  0  1  UNITS  INCH,  [hours] 

T-SCALE  =  5. 00E  +  02  UNITS  INCH,  [ft] 

KWSTflS 

RUN  2  R-TRUE  VS  TIME 


Figure  149.  Case  VI.  True  Position  Versus  Time. 
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X-SCHt_E=  1 . 00E-01  UNITS  INCH,  [hours] 
T-SCALE=5. 00E+02  UNITS  INCH,  [ft] 


KWSTRS 

RUN  2  R-HflT  VS  TIME 


Figure  ISO.  Case  VI.  Predicted  Position  Versus  Time. 
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000  002  004  006  008 


X-SCALE=  1 . OCE  +  O  1  UNITS  INCH,  [hours] 
T-SCALE=2.Q0E+01  UNITS  INCH,  [ft/sec] 

KWSTflS 

I  RUN  2  V-HRT  VS  TIME 

i. 

i 


Figure  151.  Case  VI.  Predicted  Velocity  Versus  Time. 


1 


! 


X-SCOLE=  1 .  OOE  +  0 1  UNITS  INCH,  [hours] 
Y-SCALE=2 . 00E  +  0  1  UNITS  INCH,  [ft] 


KHSTFIS 

RUN  3  R-INS  VS  TIME 


Figure  152.  Case  VI.  I.N.S.  Indicated  Position  Versus  Time. 
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X-SCALE=  1  .  OOE  +  O  1  UNITS  INCH,  [hours] 
T-SCflLE  =  5.00E  +  02  UNITS  INCH,  [ft] 


KWSTflS 

RUN  3  R-RflDRR  VS  TIME 


Figure  153.  Case  VI.  Radar  Indicated  Position  Versus  Time. 
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X-SCHLE=  1 .  OOE  +  O  1  UNITS  INCH,  [hours] 
T-SCALE=5.00E+00  UNITS  INCH,  [ft/sec] 


KWSTflS 

RUN  3  V-INS  VS  TIME 


Figure  154.  Case  VI.  I.N.S.  Indicated  Velocity  Versus  Time. 


APPENDIX  B 


SATELLITE  GEOMETRY,  VIEW  AND  RANGE  ERRORS 

The  range  measurement  equation  will  first  be  developed. 
Next  a  simple  program  has  been  formulated  to  verify  the 
"observability”  and  "suitability"  of  at  least  four 
satellites  at  any  given  time. 

A.  RANGE  MEASUREMENT  EQUATION 

The  range  measuring  process  is  characterized  by  a  set  of 
equations,  called  the  range  divergence  equations,  whicn  are 
generated  by  the  user  from  a  combination  of  I.N.S.  and 
satellite  information.  This  range  measurement  process 
involves  the  comparison  of  a  measured  value  of  range  against 
a  predicted  value  of  range. 

The  measured  range  to  a  satellite  is  determined  by 
measuring  the  incremental  phase  shift  between  the 
satellite's  clock  and  that  of  the  control  station  wnich 
supports  the  user.  These  clocks  were  synchronized  at  an 
earlier  time.  The  computed  range  on  the  other  hand  is 
obtained  from  satellite  ephemeris  data  and  user  I.N.S. 
supplied  position  information. 

The  fact  is  that  both  the  measured  and  the  computed 
range  values  contain  in  general  errors;  so  by  subtracting 
the  computed  range  value  from  the  measured  one,  the 


different  will  contain  only  the  associated  errors.  This 
difference  is  called  "the  range  divergence."  A  Kalman 
filter  can  be  constructed  to  estimate  the  errors  and  improve 
the  accuracy  of  the  raw  range  data  if  these  errors  can  be 
modeled  as  the  outputs  of  linear  systems  driven  by  white 
Gaussian  noise  [Ref.  8]* 

1 .  Range  Divergence 

The  case  of  a  single  satellite  will  be  considered, 
in  order  to  avoid  the  rotational  inconvenience  of  using 
superscripts  or  subscripts  to  keep  track  of  which  satellite 
is  being  referred  to.  The  results  are  identical  for  any  one 
of  the  18  satellites  and  therefore  very  easily  extended. 


i 


r  s  User-Satellite  position  vector 
£a  =  Earth-User  position  vector 
r_  —  Earth-Satellite  position  vector 

“ *3 

Figure  15S.  Range  Vector  Definition 
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The  range  vector  of  interest  is  the  vector  r_  from  the  user 
to  the  satellite.  It  is  explicitly  related  to  the  two 
vectors  r  and  r  which  are  defined  and  illustrated  in  the 
above  figure.  From  the  geometry  it  follows  that 

I  =  Is  -  la  (B-1> 

where 

r  =  |r|  *  |rs-ra|  =  r*r  =  (rs*ra)  •  (rg-ra)  (B-2) 

The  measured  range  to  the  satellite,  r',  is  composed  of  two 
parts 


r'  =  r  +  4r '  ( B-3  ) 

where,  r  is  the  true  range  and  Sr '  is  the  error  in  the 
measured  range  to  the  satellite.  The  computed  range  to  the 
satellite,  r",  is  in  a  similar  way  written: 

r"  =  r  +  Sr"  (B-4) 

where,  r  is  again  the  true  range  and  Sr"  is  the  error  in  the 
value  of  the  computed  range. 

The  quantity  then,  that  is  being  observed,  is  the 
difference  of  these  two  range  values  and  it  is  the  called 
"range  divergence,"  &r . 

at  =  r'  -  r"  =  Sr'  -  Sr"  (B-5) 


< 
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2.  Errors  in  Measured  Range 

To  model  the  range  measurement  error,  a  knowledge  of 
the  various  error  sources  which  are  contained  in  the 
measurement  is  required  and  fitting  these  error  sources  with 
empirical  data.  The  model  used  in  our  study  is  a  simplified 
version  of  the  one  found  in  [Ref.  4]  with  additional 
information  of  [Ref.  5].  It  is  a  linear  combination  of 
three  components  for  each  satellite  measurement  corrupted  by 
white  Gaussian  noise  (w).  Each  of  the  separate  components 
is  a  linear  system  which  is  driven  by  white  Gaussian  noise. 
The  range  measurement  error  is  modeled  by: 

ar'  =  ab  +  c ( aT  -  aT  )  +  w  (B-6) 

where 

ab  =  range  bias 
c  =  the  speed  of  the  light 
aTu  =  user  clock  phase  error 
aT3  s  satellite  clock  phase  error 
w  s  measurement  noise 


The  error  in  the  range  measurement  due  to  ionosphere  delay 
is  assumed  to  be  included  in  the  satellite  clock  phase/range 
error.  This  is  a  function  of  the  elevation  angle  and  on  the 


order  of  15  feet  as  a  good  approximation.  The  bias  terra, 
5b,  in  the  above  equation  accounts  for  the  minor  effect  of 
both  speed  of  light  bias  and  tropospheric  delay  uncertain¬ 
ties  in  each  one  of  the  four  satellite  range  measurements. 

3 .  Errors  in  Computed  Range 


The  computed  satellite  position  includes  error  which 


depends 

on  the  ephemeris  data  errors,  while  the 

I.N.3. 

errors 

account  for 

the  uncertainty 

in  the  user ' 

s  position. 

So 

far 

we  have 

r  ”  =  r  + 
-s  -s 

5£s" 

(B-7) 

r  "  s  r  + 
—a  -a 

>La" 

(3-8) 

The 

error  equation 

is  obtained  now 

since  , 

(r")2  s  r" 

.  pH 

and 

by 

taking  the 

differential  of 

both  sides  we 

get 

2r" 

5  r  "  =  r  "  •  6  r 

"  +  5  r  "  •  r  '• 

or 

5r" 

-  — 1  ( f  i*  • 

r"  - 

«r")  =  (pr  r")  * 

5r  " 

( B-9 ) 

We  can  easily  notice  that  the  quantity,  prr  ,  of  the  right 
hand  side  is  a  unit  vector  from  the  user  to  the  satellite. 

i_„  i  4r"  ;  ir„  =  [i  ,  ir  ,  i_  ]T  (B-10) 

-r  r  -  -r  rN  r£  rD 

( 
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and  also  that 


5r"  =  &r"  -  ir  "  (B-1 1 ) 

— •  — S  —a 

Finally  the  computed  error  is  written  as 

«r"  =  _ir„  •  (  5rs”  -  Sra")  (B-12) 

Without  any  significant  accuracy  loss,  it  is  assumed  that 
the  earth-satellite  range  error,  Sr_  is  approximately  zero 
for  the  following  logic  sequence.  Satellite  orbital 
parameters  are  updated  by  the  ground  tracking  network  on  a 
periodic  basis  and  relayed  to  the  user  along  with  the  range 
data.  This  ephemeris  data  is  quite  accurate  and  any 
uncer ta inties  in  computed  satellite  range  can  be  accounted 
for  by  increasing  the  satellite  clock  phase  error  [Ref.  4]. 
Thus  the  computed  range  error  can  be  written 

«r"  *  -  i^.,,  •  5£a"  (B-13) 

The  computation  of  the  above  equation  requires  values  for 
the  unit  vector  from  the  user  to  the  satellite  and  also 
current  values  for  the  north,  east  and  azimuth  I.N.S. 
position  error  states 

«£a"  =  CaN,  aE,  aD]T 
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Since  we  are  dealing  with  a  stochastic  process  simulation 
the  root-mean-sq.uared  (RMS)  values  of  the  covariance  of  the 
three  position  errors  are  used. 

The  final  form  of  the  range  divergence  equation  is 
obtained  now  by  substitution  of  equations  (B-6)  and  (B-13) 
into  the  general  form  equation  (B-5) 

Ar  =  ±r"  *  5— a”  *  26Tu  '  26Ts  *  sb  +  (3-14; 

Since  at  least  four  satellites  are  required  as  observables 
to  correct  for  the  three  components  of  position  and  the 
clock  phase  (or  time  difference),  a  minimum  of  four  range 
divergence  equations  need  to  be  solved  simultaneously. 

B.  SATELLITE  OBSERVABILITY 

m  the  following  development  of  equations  we  will  drop 
out  the  double  prime  (")  for  notational  convenience.  A 
given  satellite  must  be  in-view  by  the  user  in  order  to 
obtain  certain  measurements.  It  is  then  required  that  the 
satellite  must  be  above  some  specified  minimum  angle  of  the 
user’s  horizon  to  get  a  useful  signal.  This  minimum  angle 
depends  upon  the  capabilities  of  the  user-equipment  and  can 
be  characterized  as  arbitrary.  In  our  study  the  nominal 
value  of  this  angle  wa3  selected  of  ten  degrees.  This 
observability  criterion  together  with  the  suitable  selection 
of  18  satellites  as  the  total  number  of  satellites  for 
global  coverage,  insures  that  regardless  of  the  user’s 
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position,  a  sufficient  and  reasonable  number  of  satellites 
will  always  be  in-view,  from  which  a  "best  set”  of  the 
required  four  satellites  may  be  chosen. 

In  the  following  a  method  for  determining  whether  or  not 
a  satellite  is  observable  is  presented. 


i  satellite 


Center 

Figure  156.  Observability-Criterion  Geometry 


From  the  above  figure  and  for  the  following  calculations 
Euin  =  minimum  angle  of  elevation  for  useful  signal 


D  =  90  -  E  • 

max  7  mm 


_r  =  User-Satellite  position  vector 


r  =  Earth-Satellite  position  vector 


-3 


.1  j 

t 
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h  s  altitude  of  user 

and  the  superscript,  n,  denotes  the  frame  in  which  tne 
vector  is  coordinated  (navigation  frame). 

Since  the  earth-satellite  position  vector  in  the  earth 

•a 

frame,  r.g",  is  derived  from  the  ground  track  latitude  and 
longitude  of  the  satellite  in  the  orbit  generator  and  is 
readily  available,  the  vector  £  from  the  user  to  the  satel 
lite  coordinatized  in  the  navigation  frame  is  written  as 

0 
0 
R 
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:  ? 


The  unit  vector  along  £  is  given  by 


i  -  i  rn  • 
-r  ‘  r  L  » 


C  B— 1 7 ) 


From  the  geometry  we  observe  that  the  azimuth  component  of 
this  unit  vector  is  evaluated  as  the  sine  of  the  elevation 
angle,  E,  or  the  cosine  of  its  complementary  angle  A.  That 
is 


Ir 


r  -i 

i  0 

0 

1 


=  i  =  cos  A 
rD 


(B-18) 


3o  far  the  observability  criterion,  if  the  unit  vectors  from 
the  user  to  the  satellite  expressed  in  the  navigation  frame 
are  computed ,  becomes 


A  <  A 


or 


or 


max 


cos  A  >  cos  A 


max 


( B-1 9) 


if  ir  >  cos  >  the  satellite  is  observable 

if  if.^  £  cos  A;nax,  the  satellite  is  not  observable 

Since  we  arbitrarily  selected  for  our  study  the  minimum 
elevation  angle  of  ten  degrees  the  above  criterion  requires 
the  azimuth  component  of  this  unit  vector  to  be  greater  than 
the  value  of  co3  80°  =  o - 174- 
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ip  >  0.  174  ;  A  =  80°  (B-20) 

Total  deployment  consists  of  three  rings  or  "constella¬ 
tions"  of  six  satellites  each.  The  satellite  orbits  are 
assumed  to  be  planar  and  circular;  in  fact,  the  orbital 
speed  and  altitude  and  thus  the  orbital  period  of  all  satel¬ 
lites  is  assumed  constant.  Furthermore,  since  global  cover¬ 
age  is  desired,  the  satellites  on  any  of  the  three  rings  are 
equally  spaced;  thus,  the  circular  arc  between  any  two  adja¬ 
cent  satellites  on  a  ring  subtends  a  central  angle  of  sixty 
degrees  (6  x  60°  =  360°)  •  The  satellite  identification  code 
used  is  a  two-digit  code,  the  first  digit  denoting  the 
particular  constellation-ring  (1,  2  or  3)  and  the  second 
digit  indicating  the  particular  satellite  (1  through  6) 
among  the  six  on  the  denoted  ring.  As  an  example,  satellite 
32  is  the  second  satellite  on  the  third  ring. 

In  order  to  specify  the  orientation  of  any  one  satellite 
with  respect  to  the  Earth-fixed  frame,  three  parameters  are 
required .  The  most  convenient  parameters  are  the  Fuler 
angle,  from  which  the  direction  cosines  or  unit  vectors  to 
the  satellite  may  be  determined.  This  process  is  explicitly 
described  in  [Ref.  11]  and  here  we  will  use  directly  the 
result  for  the  unit  vector  from  Earth  to  satellite  in  Earth 
coordinates . 

-s  =  ^ a 1 3 »  a23’  a33^  (B-21) 
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where 


=  sin  c  sin  5 

a23  s  -(sin  n  cos  5  +  cos  n  sin  c  cos  K)  (B-22) 

5  con  n  cos  ?  -  sin  n  sin  c  cos  s 

and  the  three  Euler  angles  are  5,  n,  ?. 

The  initial  conditions  are  the  constant  orbital  para¬ 
meters  of  the  satellite  orbits  are  given  in  the  following 
two  tables.  Note  that  ra  refers  to  the  mth  satellite  on  the 
designated  n  ring;  for  example,  the  entry  of  2m  represents 
the  remaining  five  satellites  of  the  second  ring.  The 
missing  entries  from  the  table,  (--),  are  dependent  upon  the 
initial  values  of  the  Euler  angles  5,  n  and  c  wnicn  are 
specified.  Finally  the  latitude  and  longitude  values  refer 
to  the  ground  track  of  the  satellites. 


table  IX 

ORBITAL  DESIGN  CONSTANTS 


Orbital  Period 
Angle  of  inclination 
Altitude 

Separation  arc-angle 
Ring  spacing  arc-angle 


12  hours 

63°  (all  three  rings) 
11,000  n.  miles 
60°  (6  x  60°  =  360°) 
120°  (3  x  120°  =  360°) 
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TABLE  X 


SATELLITE  INITIAL  CONDITIONS 


Satellite  #  Initi-al  Latitude  Initial  Longitude 


1m 

21  0°  120°  (E) 

2m 

31  0°  -120°  (W) 

3m 


Angle  co 
0° 


-60(m-1 ) 

0° 

-60(m- 1 ) 
0° 


-60(m-1 ) 


Now  the  sequence  of  equations  required  to  compute  the  unit 
vectors  for  each  of  the  18  satellites  is  presented. 

5  =  63° 

n  =  1 20°  (  n-1  )  -  U)iat  (3-23 

;  =  -60(m-1 )  +  Ct 

where 


m 

n 

“ie 

C 


satellite  designator  (1  through  6) 
ring  designator  (1,  2  or  3) 

rotational  speed  of  earth  =  15°/h  =  1°/240  sec 
orbital  speed  of  the  satellite 
e  2rVn 

C  =  Ge//rs  =  — T —  =  2.  1  n. miles/sec 
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Using  the  Euler  angles  the  components  of  the  unit 
vectors  of  the  satellite-earth  range,  r se ,  are  computed 
[Ref.  10,  11]. 

a^  sin  c  sin  5 

a23  =  -(sin  n  cos  ;  +  cos  n  sin  c  cos  5)  (B-24) 

a^  cos  n  cos  c  -  sin  n  sin  c  cos  5 

The  ground  track  latitudes  and  longitudes  are  given  by 

2  2 

Latitude  8  =  arc-tan  ^a23^  +  ^a33^  ^  (B-25) 

Longitude  a  =  arc-tan  ("a23'/a33^  (B-26) 

The  required  components  of  the  unit  vector  of  the  user- 
satellite  range,  r_n ,  in  navigation  coordinates  may  no  n  be 
computed  . 


0 

'o’ 

"  r  N 

rn 

=  n3"  - 

0 

=  c  A  •  - 

e  — s 

0 

s 

rE 

R  . 

R 

r  n 

L  ■« 

L  J 

u 

-  V 

r  =  (rN)2  +  (r£)2  +  <rD)2 


(B-27) 


( B-28 ) 


ir  =  rn/r  =  [i  ,  i_  ,  ir  ]T 
“rn  “  rN  rE  rD 


(B-29) 


The  observability-criterion  requires 


5 


i  >  cos  80°  =  0. 174  (B-30) 

rD 

A  sample  output  from  the  computer  program  for  the  observa¬ 
bility-criterion  is  included  in  the  next  pages  listing  the 
satellites  which  are  in-view  at  a  particular  time  instant. 
The  output  table  presents  in  seven  columns  the  most 
important  calculated  data.  Column  1  contains  the  satellite 
number  according  to  the  previously  specified  code.  Columns 
2,  3  and  4  contain  the  north,  east,  and  azimuth  components 
of  the  unit  vector  under  consideration.  Columns  5  and  6 
give  the  ground  track  latitude  and  longitude  for  each 
satellite.  Finally  column  7  gives  the  observability  result 
denoting  with  "O.K."  the  satellites  which  fulfill  the 
criterion  and  with  " — "  those  satellites  wnich  do  not 
fulfill  the  criterion. 
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APPENDIX  D 

RICATI  COMPUTER  PROGRAM  DATA 
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